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I .  Introduction 


Background 

Research  on  efficient  negative  ion  sources  for  particle 
beam  weapons  and  advances  in  plasma  processing  of 
semiconductor  materials  have  stimulated  new  interest  in  the 
investigation  of  discharges  in  electronegative  gases.  Of 
particular  interest  to  the  Air  Force  is  development  of  high 
current  density,  low  emittance  ion  sources.  However  before 
these  plasmas  can  be  effectively  used  as  negative  ion 
sources,  the  spatial  profile  of  the  negative  ions  in 
relation  to  the  other  charged  particles  must  be  known.  A 
knowledge  of  the  spatial  profiles  will  not  only  make  it 
possible  to  efficiently  extract  negative  ions  from  the 
plasma,  but  also  allow  the  optimization  of  design 
parameters,  such  as  gas  pressure  and  temperature. 

Radio-frequency  gas  discharges  are  currently  used  to 
etch  submi rcometer  circuits  in  semiconductors  (7:77). 
Energetic  ions  that  impinge  on  the  semiconductor  material 
deliver  an  activation  energy  to  the  surface  of  the 
semiconductor  (7:77),  thereby  etching  the  material.  The 
quality  of  the  etch  is  dependent  on  the  character  of  the 
plasma  sheath,  and  the  type  and  density  of  ions  in  the 
plasma  (7:77).  In  order  that  reliability  and  consistency  be 
obtained  in  plasma  processing,  a  detailed  knowledge  of  the 
spatial  profiles  and  plasma  processes  (ion  generation/loss, 
particle  fluxes,  etc.)  is  required  (12:103). 
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Though  we  have  a  definite  need  to  fully  understand  the 
processes  in  plasmas,  our  current  knowledge  is  incomplete 
(6:145).  The  study  of  discharges  in  gases  began  as  early  as 
1924  with  the  work  of  Schottky  and  Langmuir  (2:1).  This 
early  work  concerned  discharges  in  two  component  plasmas, 
i.e.,  plasmas  comprised  of  electrons  and  positive  ions. 
Schottky's  treatment  describes  the  motion  of  charged 
particles  in  the  plasma  by  assuming  the  particle  fluxes 
result  form  a  balance  between  a  diffusive  (density  gradient) 
flux  and  a  electric  field  driven  flux.  The  predictions  made 
by  Schottky's  theory  agree  well  with  observations  (8:4699) 
Several  attempts  have  been  made  to  calculate  the  spatial 
profiles  in  three  component  plasmas  (see  References  3,8,16). 
Many  of  these  attempts  (e.g.  Reference  8)  are  extensions  of 
Schottky’s  theory  to  three  component  plasmas  (3:376). 
However,  the  results  of  these  extensions  and  other  attempts 
are  contradictory  or  of  limited  applicability. 

Problem  and  Scope 

The  purpose  of  this  study  is  to  develop  a  generalized 
model  of  a  three  component  plasma  which  predicts  spatial 
density  profiles  of  the  charged  particles.  Throughout  this 
report  the  term  model  refers  both  to  the  computational 
method  employed  and  the  relevant  plasma  kinetic  processes 


included . 


The  planned  approach  of  this  research  was  threefold. 

The  first  step  was  to  review  current  theoretical  efforts  on 
the  determination  of  the  spatial  variations  of  the  charged 
species  in  low  and  medium  pressure  gas  discharges.  This 
review  would  form  a  basis  for  determining  the  plasma 
processes  which  must  be  considered  (modeled)  and  also  guide 
the  development  of  a  computational  method.  A  further 
benefit  of  the  review  is  that  it  will  provide  sources  for 
comparison  with  theory  and  experiment  leading  to  model 
validation.  The  second  step  was  to  develop  a  model  and 
validate  it  for  plasmas  consisting  of  electrons  and  positive 
ions  only.  This  initial  model  was  to  be  developed  in  slab 
geometry  (one  dimensional)  and  later  cast  into  cylindrical 
geometry.  With  the  model  validated  for  the  two  component 
plasma  the  last  step  would  be  to  validate  the  model  for  a 
three  component  plasma. 


II.  Plasma  Theory 


In  this  chapter  we  shall  look  at  the  general  motion  of 
charged  particles  in  two  and  three  component  plasmas.  We 
will  also  look  at  the  equations  which  govern  this  motion  and 
the  relevant  production  and  loss  processes  which  occur  in  a 
plasma.  We  shall  begin  by  looking  at  the  general  physics  of 
a  plasma.  To  do  this  we  will  first  consider  a  plasma 
consisting  of  only  electrons  and  positive  ions.  After  this 
we  will  expand  the  discussion  to  a  three  component  plasma. 

The  electrons  in  a  plasma,  because  they  are  more  mobile 
(i.e.  because  of  their  smaller  mass  and  higher  temperature, 
they  have  a  larger  velocity),  will  initially  diffuse  to  the 
walls  of  a  discharge  tube  (where  recombination  occurs) 
faster  than  the  ions  (positive  ions  in  this  case). 

Therefore  the  electrons  in  the  plasma  are  lost  at  a  greater 
rate  than  the  positive  ions.  This  causes  a  charge  imbalance 
to  occur  in  the  plasma,  which  in  turn  creates  an  electric 
field.  The  field  that  develops  acts  to  reduce  the  flux  of 
the  electrons  and  increase  that  of  the  ions.  This  field 
increases  until  the  electron  and  ion  fluxes  are  equal.  When 
this  occurs,  the  sel f-cons i stant  charge  imbalance  does  not 
change  in  time  and  a  quas i - s ta t i ona ry  state  (4:197)  is 
achieved.  This  situation  is  called  arobipolar  diffusion. 

The  presence  of  negative  ions  in  the  plasma  complicates 
the  diffusion  process.  As  in  the  two  component  case,  the 
electrons  will  be  initially  lost  at  a  greater  rate  than  the 
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ions  (both  positive  and  negative).  However,  the  presence  of 
the  negative  ions  will  affect  the  net  space  charge  and  thus 
the  electric  field.  The  change  in  the  electric  field  will 


in  turn  affect  the  density  profiles  and  fluxes  of  the 
charged  particles. 

The  effect  the  presence  of  the  negative  ions  have  on 
particle  diffusion  is  critically  dependent  on  the  detachment 
rate  (15:1377).  If  the  detachment  rate  is  small,  so  that 
the  dominant  loss  of  negative  ions  is  recombination  on  the 
wall,  the  flux  of  all  three  species  is  directed  to  the  wall. 
If,  on  the  other  hand,  the  dominant  loss  of  negative  ions  is 
due  to  detachment,  the  negative  ions  have  very  little  effect 
on  the  diffusion  of  the  electrons.  In  this  case  the 
diffusion  and  density  profiles  resemble  that  of  the 
ambipolar  case  (15:1377). 

Loss  and  Production  Processes 

In  this  study  we  shall  limit  the  production  and  loss 
processes  in  the  plasma  to  four  processes.  The  first  three 
processes,  ionization,  (z),  dissociative  attachment,  (o), 
and  electron  impact  detachment,  (d),  occur  throughout  the 
discharge  volume.  Writing  the  processes  symbolically: 

e  ♦  A  -  A4  ♦  2e 

e  +  A - A* 

A*  +  e  -  A  ♦  2e 


Ionization  ( z ) 
Attachment  (a) 
Detachment  (d) 


* 


Note  that  in  defining  the  ionization  and  attachment  rates  we 
have  included  the  background  gas  (neutral  particle)  density 
in  the  expressions  z  and  a  ;  that  is,  z  =  N<?v>,  where  N  is 
the  neutral  particle  density  and  <ov>  is  the  reaction  cross 
section  for  ionizing  collisions  between  electrons  and 
neutrals.  Similarly,  <a=  N<<?v>,  where  <«?v>  is  the  reaction 
cross  section  for  attachment. 

The  fourth  process  is  recombination  at  the  discharge 
tube  wall.  These  four  processes  are  summarized  in  Table  1. 


PROCESS 

— 

Electron 

♦  Ion 

-  Ion 

Ionization 

Gain 

Gain 

Attachment 

Loss 

- 

Gain 

Detachment 

Gain 

- 

L08  8 

Recombination 

Loss 

Loss 

Loss 

at  Wall 

TABLE  1.  Summary  of  Plasma  Production  and  Loss  Processes 


Writing  the  volume  processes  in  equation  form  gives: 


(3nc/St  )  c  =  (z-o)ni  +  dnm- 
(3n-/3t)c  =  an i  -  dnin- 

(  3n  ♦  /  St )  e  -  zni 


where  ni  is 
n .  is 
n  4  is 


the  electron 
the  negative 
the  positive 


density; 
ion  density; 
ion  density. 


(  1  ) 
(2) 
(3) 
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These  equations  represent  the  local  time  rate  of  change  of 
charged  particle  density  due  to  collisions  (as  denoted  by 
the  subscript  c). 


Plasma  Equations 

The  transport  of  charged  particles  in  a  plasma  is 
governed  by  the  Boltzmann  equation  (1:206) 


3f/3t  +  v  ?f  +  (E/m)  3f/3v  =  (3f/3t)c 


(4) 


where  f(x.,v,t)  is  the  particle  distribution  function.  If  we 
assume  that  no  external  forces  act  upon  the  plasma,  the 
expression  for  the  force  (E)  can  be  written  as  E=qE»  where  E 
is  the  electric  field  generated  by  the  net  space  charge  and 
q  is  the  charge  on  the  particle  (e.g.  qs-e  for  an  electron). 
Substituting  this  into  equation  (4)  gives: 


3f/3t 


+  v  vf  +  (qE/m)3f/9v 


( 3f/3t ) c 


(5) 


Taking  the  zeroth  moment  of  equation  (5)  yields  the 
continuity  equation  (5:156): 


3n/3t  ♦  v  n<y>  =  Scon 


(6) 


where  S c o i i= ( 3n/ 3t ) c  and  <y>  is  the  average  velocity  of  the 
particle.  Scon,  therefore  is  the  density  change  due  to 
collisions,  i.e.  equations  (  1  )  —  <  3  )  .  Substituting  Scon,  for 
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each  of  the  three  species  yields  the  three  continuity 
equations : 


3nt/3t  = 

(z-a)m  +  dnin.  - 

(7  ) 

3n  ♦  /  3t  = 

zn  e  -  p  •  (^ 

(8) 

3n-/3t  = 

on  t  -  dn  in  -  -  v  •  r_ 

(9  ) 

where  £ 

is  flux  of  species  "S'1. 

The 

first  moment  of  equation 

(5)  gives  the  momentum 

equation 

(5:162): 

mn9<v>/3t  +  mn(<v>«u)<v>  +  kTvn  -  q£n  =  mn^T  (<v'l>-<v>) 

-  m<v>Scoi.i  (  10) 

where : 

v'  is  the  momentum  transfer  frequency  (collisions 

H  T 

with  neutrals ) ; 

<vN>  is  the  average  velocity  of  neutrals,  which  is 
assumed  to  be  zero  (the  neutral  particle  velocity  is  assumed 
Maxwel 1 ian ) ; 

k  is  the  Boltzmann  constant  and  T  is  the  species 
temperature.  In  deriving  this  expression  for  the  momentum 
equation,  it  is  assumed  that  the  pressure  tensor  can  be 
represented  as  a  scalar  pressure  (kTvn). 


Equations  (7)-(9)  and  (10)  form  the  basic  equations 
needed  to  model  plasma  behavior.  These  equations  are  then 
complemented  with  Poisson’s  equation: 


v£  =  (e/£#Kn*-ni-n.) 


(11) 


III.  Review  of  Plasma  Models 


Several  attempts  to  model  plasma  behavior  have  been 
made.  In  this  section  we  shall  review  four  models  in  an 
attempt  to  identify  critical  assumptions  and  demonstrate  the 
limitations  and  the  disparity  in  the  results  of  these 
approaches.  The  four  models  which  shall  be  reviewed  are: 

(1)  Ambipolar;  (2)  Self  and  Ewald;  (3)  Lee  and  Lewis;  (4) 

Von  Engel  and  Edgley. 


Ambipolar  Model 

The  ambipolar  model  is  the  simplest  of  the  models.  In 
this  model  it  is  assumed  that  the  plasma  has  reached  a 
steady  state;  therefore,  the  time  dependence  of  the 
continuity  and  momentum  equations  is  eliminated.  The 
momentum  equation  for  the  electrons  becomes  (in  one 
dimension ) : 


n  tm<v>  3<v>/dx  ♦  kTjSm/; 


(  12  ) 


Substituting  for  the  particle  charge  (q=-e)  and  assuming 
that  the  inertial  term  (<v>3<v>/9x)  can  be  ignored  due  to 
the  small  mass  of  the  electron: 


kT*3ni/3x  +  nieE  =  -nm^  <v> 
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(13) 


This  equation  can  be  used  to  solve  for  the  ambipolar 
drift  velocity,  <v>.  To  do  this  we  first  define  nobility  as 


*/=|q|/mv  (1:137).  We  will  also  define  a  diffusion 

coefficient  using  the  Einstein  relation,  D  =  kT/mv>  (1:137). 
Using  these  definitions  in  equation  (13),  one  obtains: 

<v>  =  -  (  D i /n  f  )  3n  i /9x  - 

But  flux,  PE  is  equal  to  nt<v>,  therefore: 

rE  s  -Di9ni/Bx  -  *  Ens  (14) 

Assuming  the  positive  ion  flux  can  be  written  in  the  Bane 
fashion,  we  have: 


^  =  -D«3n«/3x  -  M  +  En ♦  (15) 

It  is  next  assumed  that  the  charged  particle  density 
profiles  are  approximately  the  same;  therefore,  n*=n«=n 
(1:137).  In  the  steady  state,  the  electron  and  ion  fluxes 
are  equal.  Therefore,  by  setting  equations  (14)  and  (15) 
equal  to  each  other  and  assuming  quas i -neutral i ty  (n«=n,=n) 
it  is  possible  to  derive  an  expression  for  the  electric 
field.  The  resulting  expression  (1:138): 


E  =  { (  D  ♦  -  D  i )  /  ( >u+  +  </E)J  vn/n 


(  16  ) 


Substituting  equation  (16)  into  equation  (14)  results  in 
(1:138): 


T  =  -D*5n/3x 


(  17  ) 


where  Da  is  the  ambipolar  diffusion  coefficient; 

D*s(i/tDi  +  *[D*)/(«4+  u  £  )  . 

Substitution  of  equation  (17)  into  equation  (7), 
(noting  that  3n/3t  =  0  in  steady  state),  we  obtain  ( a  and  d 
are  zero  for  a  two  component  plasma): 


DA3?n/3x*  +  zn  =  0 


(  18) 


Boundary  conditions  must  be  established  before  this 
differential  equation  can  be  solved.  We  expect  that  the 
discharge  is  symmetric  about  the  center  of  the  discharge 
tube  (two  parallel  plates  separated  by  a  distance  2L  in  our 
case);  therefore  the  flux  a  the  center  (x=0)  should  be  zero. 
If  the  flux  is  zero  at  the  center  of  the  discharge,  then  by 
equation  (17)  ( 3n/5x ) * r o=0 .  Ignoring  the  development  of  a 
sheath  area,  it  is  assumed  that  the  particle  density  at  the 
discharge  wall  (x=L)  is  zero.  Applying  these  boundary 
conditions,  equation  (18)  is  easily  solved  to  obtain  (for 
rectangular  1-D) : 


n  =  n c oCos ( *x/2L ) 


(  19) 
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where  l z/Da ] 1  ,=*/2L. 

Substituting  equation  (19)  into  equation  (16)  gives: 

E  =  (  ( D  ♦  -  D«)/Uf  ♦  )  )  (  2L/*)Tan(tix/2L)  (20) 


An  example  of  the  ambipolar  density  profile  is  given  in 
Figure  1.  Though  in  our  development  we  assumed  the  particle 
profiles  to  be  the  same,  the  figure  illustrates  the 
difference  in  the  two  profiles  which  must  exist  to  develop 
the  field  given  by  equation  (20). 


Figure  1  Charged  Particle  Profiles  Under  Ambipolar 
Conditions  (Golant:198) 
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the  ambipolar  model. 


The  first  difference  between  this  model  and  the 
ambipolar  model  can  be  seen  in  the  momentum  equations.  Self 
and  Ewald  make  the  same  assumptions  in  developing  their 
electron  momentum  equation  as  are  made  in  the  ambipolar 
model,  but  unlike  the  ambipolar  model,  they  do  not  assume 
the  electron  equation  can  also  be  applied  to  the  positive 
ions.  In  their  development,  Self  and  Ewald  retain  the 
inertial  term  (<v>9<v>/9x)  in  the  ion  momentum  equation. 

An  additional  important  difference  in  Self  and  Ewald’s 
model  is  that  they  do  not  assume  a  zero  density  at  the 
discharge  wall.  They  instead  attempt  to  model  the 
plasma-sheath  boundary  by  applying  the  Bohm  criterion 
(14:2487) . 

With  the  major  differences  in  this  model  and  the 
ambipolar  model  identified,  an  outline  of  Self  and  Ewald’s 
model  can  be  discussed  (a  detailed  development  of  Self  and 
Ewald's  model  is  included  as  Appendix  A).  Like  the 
Ambipolar  model,  Self  and  Ewald  assume  quasi-neutrality 
(nc=n,=n).  This  reduces  the  problem  to  one  of  solving  three 
equations,  a  continuity  equation  and  two  momentum  equations. 
Normalization  of  the  continuity  and  momentum  equations 
results  in  equations  ( A- 1 2 ) - ( A- 1 4 )  (in  one  dimension). 

d(su)/ds  ud{Ln(n)]/ds  =  A  (  A  —  1  2  ) 

u  =  R’(d'i/ds)-  d[Ln(n))/ds  ( A-  13) 

u  +  u(du/ds)  =  dn/ds  -  x  d[Ln(n)]/ds  < A  —  1 4 ) 
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where : 


A  is  the  ratio  of  the  ionization  rate  to  the  sum  of  the 
momentum  transfer  collision  frequency  and  ionization  rate; 
u  is  the  normalized  particle  velocity; 
s  is  the  normalized  spatial  coordinate; 
n  is  the  normalized  electric  potential; 

I  is  ratio  of  the  ion  temperature  to  the  electron 
temperature ; 

R’  is  a  normalizing  factor  ( m «  c+ /m  e  Og )  • 

Making  the  normalized  particle  velocity  the  independent 
variable,  Self  and  Ewald  solve  for  ds/du,  s,  Ln(n/no),  and 
the  normalized  potential  (n)  (no  is  the  particle  density  at 
x=0).  The  resulting  equations  are: 


ds/du  =  { ( l+T)-u2)/[A( l+t)+  R'u2]  ( A- 20 ) 

s=  [  (  1+1) ‘ ' * (R"+A) ] / [ A  * / * <R- ) 2 '2 ]Tan- 1 ( Wu ) -  u/R’  (A-21) 

Ln ( n/n  o  )  =  - { [ R ' +A ] / ( 2R ' ) } Ln ( Y )  (A-22) 

’t  =  -{  (  1-X/R'  )  (R'+  A  )  /  [  2  (  R  "  )  *  ]  }  Ln  (  Y  )  -u  2  /  [  2  (  1 +R  ’  )  ]  (A-23) 


Where 

R"=  ( 1  +  1/R  1  )  ; 

W  =  [ ( 1+1/R' ) / A ( 1 + I ) ) » - 2] ; 

Y  =  { (R"+A]u2+A< 1+T) }/[A( 1+T) ] . 

Equation  (A-20)  is  used  to  place  an  upper  bound  on  the 
variable  u.  The  upper  limit  of  u,  is  the  value  which  makes 
equation  (A-20)  equal  to  zero.  When  this  condition  is 


f 


1  5 


reached,  the  ion  velocity  is  equal  to  the  ion  sound  speed 
(the  Bohm  criterion)  (14:2487).  The  discharge  boundary  is 


the  value  of  "s"  when  the  upper  value  of  u  is  substituted 
into  equation  (A-21)  (14:2488).  Therefore  in  this  model  the 
density  is  not  constrained,  as  it  was  in  the  ambipolar  model 
(i.e.,  the  density  made  to  be  zero  at  the  discharge 
boundary ) . 

The  density  profile  at  the  discharge  boundary  is 
dependent  on  the  pressure  of  the  background  gas.  Self  and 
Ewald  incorporate  this  pressure  dependence  in  their 
equations  by  use  of  the  pressure  parameter  (A)  (14:2488). 

In  low  pressure  plasmas  the  ionization  rate  is  much  greater 
than  the  ion  collision  frequency  (14:2487)  (mean  free  path 
between  collisions  is  large);  therefore,  A  approaches  unity. 
In  these  cases  the  density  profiles  approach  that  of  the 
Free-Fall  theory  (14:2488).  Free-Fall  theory  assumes  that 
that  the  motion  of  the  charged  species  is  not  impeded  by 
collisions.  The  velocity  of  the  particles  in  this  theory  is 
equal  to  that  the  particles  would  have  by  falling  through 
the  electric  potential  existing  in  the  plasma.  As  A 
approaches  zero  (high  pressure  plasmas)  the  density  profile 
transitions  to  the  results  of  the  ambipolar  model  (14:2488). 
Figure  2  depicts  this  transition  from  the  free  fall  theory 
to  the  ambipolar  results  for  increasing  A. 
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Figure  2  Density  Profiles  for  Plane-Symmetric  Discharge 

(Self : 2 4 88  ) 


Lee  and  Lewis’  model  was  originally  developed  for 
discharges  in  oxygen;  however,  given  the  ionization, 
attachment,  and  detachment  rates  for  a  particular  gas,  the 
model  can  be  applied.  A  detailed  analysis  of  this  model  was 
made  in  reference  2. 

We  will  first  examine  the  assumptions  Lee  uses  to 
de ve 1  op  his  solution.  As  before,  we  shall  begin  with  the 


continuity  and  momentum  equations.  Lee  assumes  steady  state 
solutions;  therefore,  the  continuity  equations  <  7 ) - <  9 ) 
reduce  to  (in  cylindrical  geometry): 


r-*9(rrf  )/9r  -  (z-a)nt  -  dntn.  :  0  (21) 

r-'S(rr4  )/5r  -  zn  e  =  0  (22) 

r-*9(rP_  )/3r  -  ons  +  dntn-  =  0  (23) 

where  r  is  the  radial  coordinate  in  cylindrical  geometry. 

As  in  the  ambipolar  model,  Lee  assumes  the  fluxes  of 
each  species  can  be  written  as: 

fs-Dtnt/r-^ntE  (24) 

r+  =  -D.  n,/  r  +  ^n«E  (25) 

r  =  -D-  n-/  r  -  i/_n-E  (26) 

Using  that  in  the  steady  state,  P  =  +  and  Poisson’s 

equation,  Lee  is  able  to  eliminate  the  electric  field  from 
these  equations.  The  end  result  is  a  set  of  four 
differential  equations  for  the  electron  and  negative  ion 
fluxes  and  densities. 

x  *  1  ( xy  s )  =  [(z-o)/dyi  +  (d/o)yj  (27) 

x  *  1 ( xy « )  =  yi-  (d/a)yi  (28) 

y  i  •  +  2 ( T./T £ ) y 2 -  Ay  3  -  By «  (29) 

yi’y*  -  (T-/Tt )y,  ’yi=  (B-A)yty«  -  Cyjyj  (30) 
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x  is  the  normalized  radius  (r/R)  (R  is  outer  tube  radius) 
yi  is  the  normalized  electron  density  ( n t ( r ) /n t < 0  )  ) 
y*  is  the  normalized  negative  ion  density  ( n - ( r ) /n i ( 0  )  ] 
y3  is  a  normalized  electron  flux  [  (  r ) /.3n  c  (  0  )R ) 

y«  is  a  normalized  negative  ion  flux  [  (r)/oni(0)R) 

A  -  m «  O  +  aR 2 /T  t ;  B-A+(m-v_aR*/Ti);  C -  mi  joRJ/T i 

This  set  of  four  differential  equations  (eq.  (27)-(30))  is 
solved  using  a  fourth  order  Runge-Kutta  method  (2:24). 

Like  the  ambipolar  model,  Lee  ignored  the  sheath  region 
and  assumed  the  density  of  all  three  species  to  be  zero  at 
the  wall.  Lee's  second  boundary  condition  was  that  the 
particles  fluxes  at  the  axis  were  zero  because  of  symmetry. 

In  testing  this  model  the  primary  problem  was  the 
sensitivity  of  the  calculation  and  subsequent  convergence  to 
the  selected  value  of  the  on-axis  ratio  of  negative  ions  to 
electrons  (h).  This  parameter  was  the  starting  value  used 
in  the  numerical  integration.  The  convergence  to  a  solution 
was  slow  for  a  poor  choice  of  h.  This  sensitivity  to  the 
parameter  h  is  discussed  in  detail  in  reference  2. 

An  example  of  the  electron  and  negative  ion  profiles 
predicted  by  Lee’s  model  is  shown  in  Figure  3.  The  profiles 
in  Figure  3  are  representative  of  the  profiles  presented  in 
Lee’s  paper  (Ref.  8)  and  Clouse’s  analysis  (Ref.  2).  In 
particular  we  should  note  from  figure  3  that  Lee’s  code 
predicts  that  the  negative  ions  are  distributed  throughout 


This  model  resembles  Lee’s  model  in  many  ways. 


As 


before  we  shall  proceed  by  discussing  the  development  of  the 


continuity  and  momentum  equations,  and  then  discuss  the 


boundary  conditions  applied  in  the  model 


The  particle  continuity  equation  and  the  electron 


momentum  equation  used  by  Von  Engel  are  the  same  as  those 


used  by  Lee.  The  ion  momentum  equations  however  differ 


significantly.  In  developing  his  positive  ion  equation  Von 


Engel  assumes  that  the  electric  field  driven  drift  velocity 


will  be  much  greater  than  the  positive  ions’  thermal 


velocity;  therefore,  the  term  involving  the  ion  temperature 


can  be  ignored.  Making  this  assumption  leaves  us  with  the 


following  (in  cylindrical  geometry): 


n»3<v,>/5r  =  -eE/(m«<v«>)  -  v>  +  -  zni 


(  31  ) 


For  the  negative  ions  Von  Engel  assumes  that  the 


negative  ion  velocity  will  be  primarily  thermal.  He  further 


assumes  the  the  motion  of  the  negative  ions  will  be  governed 


by  "the  balance  between  of  the  outward  diffusion  due  to  the 


density  gradient  and  the  inward  drift  due  to  the  electric 


field"  (3:378).  Using  these  assumptions,  Von  Engel  drops 


the  inertial  term  from  the  negative  ion  momentum  equation 


This  reduces  the  negative  ion  momentum  equation  to: 


5n-/9r  =  en-E/(kT.)  -  (  V. m - / ( kT  -  )  ] n - < v . > 


(  32  ) 


•  mj>  ’ji  *  a  *  .  *  - 


Like  Lee,  Von  Engel  uses  the  continuity  equations, 
■omentum  equations,  and  Poisson’s  Law  to  develop  a  set  of 
differential  equations  which  he  solves  using  a  fourth  order 


% 

Runge-Kutta  method. 

The  other  difference  between  Lee’s  and  Von  Engel  models 
is  the  boundary  condition  at  the  wall.  Where  Lee  assumes  a 
zero  density,  Von  Engel  assumes  a  condition  much  like  that 
of  Self  and  Ewald’s  model.  Von  Engel  terminates  his 
numerical  integration  at  the  position  where  the  following 
boundary  condition  for  the  electrons  is  met  (3:379): 


ue  -  (vi>/vje-  ( 2m * /nm t ) * 1  * 


With  this  boundary  condition  Von  Engel  assumes  that  the 
discharge  boundary  is  the  point  where  the  normalized 
electron  velocity,  ui,  is  equal  to  the  ratio  of  the  electron 
thermal  velocity  to  the  ion  sound  speed  (3:379).  In  this 
way,  Von  Engel  is  able  to  place  a  boundary  condition  on  the 
electrons  without  constraining  the  density.  Von  Engel 
believes  that  forcing  all  particles  densities  to  be  zero  at 
the  wall  over  constrains  the  solution  (3:376). 

Von  Engel  noted  computational  problems  similar  to 
Lee’s  in  that  the  process  is  slow  and  sensitive  to  the 
choice  of  h  (3:381).  An  example  of  Von  Engel’s  results  are 
presented  in  figure  4. 
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norma  1 i zaton ,  Lee  used  n - / n  -  *  (where  n-o  is  the  negative  ion 
density  on  axis)  while  Von  Engel  used  n»-/nto-  The  most 
striking  difference  in  the  two  author's  profiles  is  that 
Von  Engel  predicts  that  the  negative  ions  are  confined  near 
the  center  of  the  discharge,  while  Lee’s  negative  ions  are 
distributed  throughout.  In  Figure  5  we  have  attempted  to 
demonstrate  this  difference  by  placing  Von  Engel’s  negative 
ion  profile  on  Lee’s  plot  (Figure  3).  This  difference  may 
be  related  to  a  difference  in  detachment  rate  (as  reported 
by  Tsendin).  Unfortunately  this  could  not  be  determined 
from  the  information  given  in  the  authors’  papers. 


IV.  Model  Development 


The  numerical  method  employed  is  based  on  a  finite 
difference  solution  of  continuity  and  momentum  equations  of 
the  three  charged  species.  One  advantage  of  this  method  is 
that  it  would  allow  a  time  development  of  the  density 
profiles  and  electric  field.  This  removes  the  sensitivity 
of  the  solution  from  the  choice  of  initial  conditions  (which 
was  seen  in  the  solutions  of  Lee  and  Von  Engel).  Because 
the  density  profiles  are  allowed  to  adjust  themselves  in 
time,  our  initial  profiles  can  evolve  into  the  steady  state 
solution,  as  opposed  to  requiring  an  iterative  selection  of 
new  initial  conditions  when  the  solution  does  not  meet  the 
boundary  conditions  in  the  steady  state  (as  in  Lee  and  Von 
Engel’s  solutions). 

The  discussion  in  this  chapter  is  divided  into  three 
parts.  The  first  of  these  describes  the  general  finite 
differencing  method  (as  applied  to  a  flux  divergent 
problem).  The  second  part  will  cover  the  development  of  the 
finite  differenced  plasma  equations  (continuity  and 
momentum).  The  discussion  in  this  chapter  ends  with  a 


general  description  of  the  model  code. 


however  the  general  principles  are  the  same  for  all 
geometries.  As  in  any  finite  differencing  scheme,  space  is 
divided  into  cells.  In  this  case,  "space"  is  the  discharge 
area  between  one  wall  of  the  discharge  tube  and  the  center 
(because  the  discharge  is  symmetric  about  the  center  only 
one  side  of  the  discharge  need  be  modeled). 

To  explain  the  setup  of  the  finite  differencing  grid 
and  general  principle  of  the  finite  differencing  method,  we 
shall  consider  a  simple  one-dimensional  flux  divergent 
problem.  For  this  example,  assume  that  "space"  is  divided 
into  five  cells  (see  Figure  6).  Each  cell  (i),  has  a 
density  m ,  defined  at  the  cell  center.  Fluxes  are  defined 
at  the  cell  boundaries. 


Looking  now  at  a  single  cell,  we  find  that  the  change 
in  the  number  of  particles  in  unit  time  is  given  by: 

{change  in  density/unit  timeJdVol  =  {Net  Flux  into  CellfdA 
♦  { Product i on/Losses ) dVol 

In  this  example,  we  shall  assume  that  there  are  no 
production  or  loss  mechanisms  within  the  volume  of  the  cell 
therefore,  the  change  in  the  number  of  particles  in  a  cell 
will  be  determined  by  the  fluxes  into  and  out  of  the  cell. 
Putting  our  conservation  statement  into  equation  form: 

{ An/At ) AXAyAZ  =  -<r  -  r  )a>az 

Cut  in 

This  is  actually  the  continuity  equation  (eq.  6)  with  the 
collisional  production  (Seoul  term  equal  to  zero; 
therefore: 

-n/~\  -  -  ~ir  )/dx 

Appr ox  1  mat  1 ng  the  time  derivative  with  a  forward  differenre 
gives: 

{  n  (  t  +  i )  -  n  (  t  )  }  / 1  r  -C0ut  -  r,n  )/dx  (33 

where  I  is  the  time  step  interval.  From  Figure  6  we  see 
that  for  the  first  cell,  the  in  flux"  is  ~  and  the  "out 


flux'*  is  .  Letting  h  =  dx  and  using  an  index  1  to  denote 
cell  numbers,  equation  (33)  becomes: 

(n,(t  +  i)  -  n,(t)  )  / 1  =  -(q  -  r,_j  )/h  (34) 

Equation  (34)  can  be  solved  for  n,(t+t)  to  give: 

nVJ=  nji  -  (l/h)(rj  -  r3  )  (35) 

where  j  denotes  the  current  time  step  (t)  and  j+1  denotes 
the  next  time  step  (t+T). 

So  far  in  our  discussion  we  have  assumed  that  the 
finite  difference  equations  (e.g.,  equation  (35))  are 
numerically  stable.  This  will  only  be  true  if  the  Courant 
condition  is  met.  The  Courant  condition  is  (10:626): 

I  v  c I  I / h  s  1  (36) 

where  vc  is  the  characteristic  speed  at  which  information 
travels  through  the  finite  difference  grid.  In  the  example 
above,  the  characteristic  speed  is  the  maximum  velocity  of 
the  particles.  The  necessity  of  obeying  the  Courant 
condition  can  be  better  understood  at  examining  the  way  in 
which  the  density  of  a  cell  is  calculated.  The  cell  density 
is,  in  part,  calculated  from  information  (fluxes)  gained 
from  the  boundaries  of  that  cell.  If  the  time  step  (x)  is 


to  large  (i.e.,  vc:  >  h  ),  the  information  (the  fluxes  in 


our  case)  required  by  the  differencing  scheme  lies  at  points 
outside  the  boundaries  of  the  cell.  Since  the  fluxes  will 
not  be  defined  at  these  points,  a  lack  of  information 
exists.  This  lack  of  information  gives  rise  to  an 
i  nstabi  lity  (10:627).  Therefore,  provided  that  (i/h)s(  1/vd 
the  Courant  condition  assures  us  that  our  finite 
differencing  scheme  will  remain  numerically  stable. 

We  now  increase  the  complexity  of  our  example  by 
including  the  influence  of  an  electric  field.  In  this  way, 
we  can  examine  the  use  of  the  finite  differencing  scheme  to 
calculate  plasma  parameters  (in  a  simple  plasma).  For  this 
example  we  shall  assume  ambipolar  conditions;  therefore,  the 
electron  continuity  equation  and  flux  are  given  by: 

3n/3t  =  - ( 3r/ 3x )  +  zn  (37) 

r'  =  -  D  5n  /  3x  -  „  E  7n  (38) 

(For  the  purpose  of  this  example  only  the  development  of  the 
electron  equations  will  be  discussed).  Finite  differencing 
of  equations  (37)  and  (38)  results  in  the  following: 

n:.+  1=  n'\  -  (T/h)(rj  -  ly  )  ♦  zin’i  (39) 

ry  =  -  (  D  i  / h  )  (  n"'i  ♦  i  -  n\  )  -  (  „  e\  /  2  )  (  nJ,  .  i  ♦  nJ,  )  (40) 

Before  proceeding,  a  few  notes  on  equation  (40)  are  in 
order.  The  fluxes  in  the  finite  differencing  grid  are 
calculated  on  the  cell  boundaries;  therefore,  we  shall 


calculate  the  density  and  gradient  terms  in  equation  (40)  at 
the  cell  boundaries.  Because  the  densities  are  defined  at 
the  cell  centers,  the  density  at  the  cell  boundary  is 
approximated  by  averaging  the  densities  of  the  cells 
adjacent  to  the  boundary  (thus  giving  the  [ni+nml/2  term). 

Before  equations  (39)  and  (40)  can  be  developed  into  a 
set  of  equations  which  can  be  used  to  calculate  the  cell 
densities,  boundary  conditions  must  be  established.  For 
this  example,  it  will  be  convenient  to  establish  boundary 
conditions  on  the  flux  at  the  discharge  axis  (x  =  0)  and  the 
discharge  wall  (x  =  L).  Assuming  that  the  discharge  will  be 
symmetric  about  the  discharge  axis,  the  flux  at  x=0  is  zero 
(i.e.  r  =  0).  For  the  wall  condition  we  shall  assume  the 
flux  at  the  wall  is  equal  to  a  constant  (i.e.  rE  =  L). 

Substituting  the  flux  expression  (40)  into  (39)  results 
in  an  expression  which  can  be  used  to  solve  for  the  electron 
density  in  cells  2  thru  N-l ,as  shown  below: 

l((Dt/hJ)(n,',-i-2ni  +nJ,  .i)-(^[/2l((nJl,1  +n'1,  )  E'l 
-  (  n', -n  ,J  -  i  )  E  ,J  -  i  )  }  ♦  zin  i  2  £  i  s  N-l  (41a 

Incorporating  the  boundary  conditions  for  rj  and  TN  in 
equation  (39),  yields: 

n  j  =  nJi  -H(Di/h*)(n?«i-  nji)/h  -  [u  E  i/2  )  (  n°I  .  i  +  nj,  )  )  ♦  izn'i 

i  =  1 


(4  1b 


n  m  = 1  n«  -  X  |  L  -(Di/h1  I  ln«-  n»  -  i  )  -  (  ^  -  i  /  2  )  (n'ii  +n«-i))+Tz  n'» 

i=N  ( 4 1 c  > 

A  similar  set  of  equations  can  be  developed  for  the  positive 
ions. 

The  electric  field  is  calculated  using  Gauss’  Law. 

t 

{  E  (  x  )  -  E  (  0  )  }  =  (  e  /  E  o  )  ^  (  n  ♦  -  nE)  dx' 

c 

Approximating  this  integral  with  a  summation,  the  electric 
field,  directed  toward  the  boundary  is: 

N 

E,  =  X  (  n  »  ,  t  -  n  i  ,  r  )  h 

r  *  1 

where  dx  has  been  replaced  by  h  and  due  to  symmetry,  Eo  =  0. 

Summarizing  the  procedure  used  in  calculating  the 
density  profiles: 

1)  Start  with  a  set  of  initial  conditions  for  the 
electric  field  and  electron  and  positive  ion  densities. 

2)  Calculate  the  electron  and  ion  densities  at  an 
advanced  time  using  the  density  equations. 

3)  Compute  the  electric  field  in  this  advanced  time 
using  the  new  densities. 

4)  Increment  time  (i.e.  redefine  j+1  values  as  j 
values),  repeat  until  process  converges. 

Once  again  the  Courant  condition  (equation  (36))  roust 
be  obeyed  for  the  differencing  scheme  to  remain  stable. 


However,  the  speed  at  which  information  travels  through  the 
grid  is  complicated  by  the  electric  field.  The  electric 
field  term  in  the  density  equations  (equations  (41a)  - 
(41c))  is  itself  a  function  of  density.  This  situation 
produces  a  feedback  which  acts  to  make  the  differencing 
scheme  more  sensitive  to  the  time  step  than  in  the  previous 
case  (simple  diffusion,  no  electric  field).  This  feedback 
phemonena  will  be  discussed  in  greater  detail  at  a  later 
time.  However,  it  should  be  possible  to  put  limits  on  the 
time  step  required.  The  speed  of  information  in  the  grid 
should  lie  between  the  maximum  velocity  of  the  particles  and 
the  speed  of  light,  c  (the  speed  at  which  changes  in  the 
electric  field  occur).  Therefore  the  upper  limit  on  the 
time  step  (the  largest  possible  time  step)  can  be  found  by 
using  the  particle  velocity  in  equation  (36).  The  lower 
limit  of  the  time  step  is  calculated  by  substituting  c  into 
equat ion  ( 36 ) . 


Devel opment  o f  Mode  1  Eauat  ions 

We  are  now  ready  to  begin  the  development  of  the 
equations  which  we  will  use  to  calculate  the  density- 
profiles  in  a  three  component  plasma.  This  discussion  of 
the  model  equations  development  is  divided  into  three  areas: 
1)  boundary  conditions;  2)  momentum  equations;  3)  difference 
equations.  For  simplicity,  we  will  develop  the  equations  in 
a  one  dimensional  (slab)  geometry;  however  the  development 
can  be  easily  applied  to  other  geometries.  Therefore  in  our 


development  we  will  assume  that  the  plasma  is  confined 
between  to  infinite  parallel  plates  a  distance  2L  apart. 

Boundary  Cond i ti ons .  As  in  the  example  above  boundary 
conditions  are  applied  to  the  particle  fluxes.  As  before, 
we  shall  assume  that  the  discharge  is  symmetric  and 
therefore  the  flux  of  all  three  species  is  zero  at  the 
center.  The  flux  cond ition  on  the  wall  for  the  electrons 
and  negative  ions  is  derived  from  basic  transport  theory 
(9:496).  According  to  transport  theory,  the  thermal  flux  on 
a  completely  absorbing  boundary  (which  the  wall  is  assumed 
to  be)  is  given  by  <l/4)n(L)<u>  (9:496).  Where  <u>  is  the 

average  speed  of  the  particle  and  n(L)  is  the  density  at  the 
boundary . 

In  deriving  the  electron  flux  boundary  condition,  it  is 
assumed  the  the  average  velocity  of  the  electrons  incident 
on  the  wall  is  given  by  their  thermal  velocity.  Ising  the 
thermal  velocity  predicted  by  kinetic  theory, 

<u) : ( 8k  I  1  1  (1  3:268),  results  in  the  following 

boundary  condition  on  the  electron  flux: 

:  s  m(L)(kTi/2H«  »  (42) 


A  phenomenon  expression  for  the  wall  flux  of  the 
negative  ions  can  be  derived.  Since  the  effect  of  the 
electric  field,  is  to  accelerate  the  negative  ions  toward 


the  center  of  the  discharge  (where  the  field  is  weaker);  the 


characteristic  velocity  governing  the  negative  ion  flux  at 


the  wall,  is  the  ions  average  thermal  velocity.  Therefore, 
we  shall  assume  drift  velocity  of  the  negative  ions  is  small 
compared  to  their  thermal  velocity.  Therefore  the  wall  flux 
for  the  negative  ions  is  given  by: 


r  =  ( n  -  v  T  h ) / 4  =  n  -  (  L  ) (  k  T  -  /  2  *  )  1  2  (43) 

In  the  case  of  the  positive  ions  the  Bohm  criterion  is 
applied.  The  Bohm  criterion  states  that  the  ion  velocity  at 
the  discharge  boundary  must  be  greater  than  or  equal  to  the 
ion  sound  speed.  The  ion  sound  speed  is  given  by  (14:2487): 

<u>=  [k(Tt  +  T, )/m, ]  1  *  (44  ) 


Because  flux  is  equal  to  n<u>,  the  following  expression  of 
the  positive  ion  flux  at  the  wall  is  obtained  (note:  we  are 
not  assuming  a  thermal  flux  and  therefore,  do  not  use 
(n<u>)/4  as  the  wall  flux  expression): 

r*  =  n*(L)(k(Ti  +  T.l/m,)1  2  (45) 


Momentum  Equations.  With  the  derivation 
boundary  conditions  out  of  the  way  we  can  now 
at  the  development  of  the  momentum  equal  inns, 
with  the  general  momentum  equation  (10): 


of  t  he 

begin  to  look 
We  will  start 


(  10) 


mn9<v>/9t  ♦  mnl  <v>»v)<v>  ♦  kTtzn  -  qEn  =  -ran.  <\ > 

-  _  MT  - 

-  ID  <  V  >  S  c  «  I  1 

We  will  first  consider  the  electron  momentum  equation. 
Like  all  the  models  reviewed,  the  inertial  and  collisional 
terms  will  be  ignored.  For  a  quasi-steady  state  the 
electron  momentum  equation  yields: 


r  -  D  e  9n  e  /  9x  -  p  En  e 
E  E 


(46) 


The  positive  ion  momentum  equation  proved  to  be  more 
complicated  than  the  electron  equation.  The  larger  mass  of 
the  ion,  coupled  with  a  significant  drift  velocity  (due  to 
the  electric  field)  makes  it  impossible  to  ignore  the 
inertial  and  collisional  terms.  Because  we  expect  the 
velocity  of  the  positive  ions  to  be  mainly  dependent  on  the 
electric  field,  the  diffusive  flux  term  (kT9n/9x)  can  be 
ignored.  Ignoring  the  diffusive  term  and  substituting  njz 
for  Scou  reduces  the  steady  state  momentum  equation  for 
positive  ions  to: 


<v,>?<v,>/5x  =  (eE/m, )  -  v  <  v ,  >  -  (ntz/n«)<v«>  (47) 


Equation  (45)  can  be  used  to  find  an  expression  for 
calculating  the  positive  ion  velocity. 


■ . 

v./. 
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The  mean  free  path  for  collisions  between  neutral  gas 
atoBB/noleculeB  and  positive  ions  is  given  by: 

A  =  1  /  <N<?  > 

where : 

N  is  the  neutral  gas  density 
o  is  the  cross  section  for  collisions. 

Using  the  definition  v=No<v>  the  expression  for  < v ♦ >  can 
be  rewritten  in  the  following  way: 

<v«>  -  No<v»>*  s  <v«>*/X  (48) 

Substituting  equation  (48)  into  equation  (47)  and  rewriting 
<v*>9<v4>/3x  as  9(<v«>*/2)/3x  yields: 

9(<v«>J/2)/9x  =  eE/nu-  <v  ♦>*/  ».  -  <v«>znt/n.  (49) 

Finite  differencing  equation  (49)  gives: 

( v?,  -  v?t-i)/(2h)  =  eE,/m.  -  v,  /  •  -  zv,(m/n«) 

where  (m/n«)=  (nc.m  +  n  i  .  ,  )  /  (  n  .  .  ,  «  i  +  n.,t).  This  is  the 
ratio  of  the  electron  to  ion  density  (note  the  densities 
were  averaged  to  approximate  the  density  at  the  cell 
boundary  were  the  flux  and  velocity  are  defined). 

Collecting  like  terms  and  using  the  quadratic  equation  to 
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solve  for  v,  results  in: 


vf=  (  -  l  -  ’  z  (  nl/ri’*  )  1  ♦  (  I  v  '  z  ( n  i/nt  )  ]  *+  4<A'(eE/n  ♦ 

vv]  .  i/2h)  >  *  /  *)/2 


(50) 


where  a’=  2h> / { A+2h ) .  Using  the  boundary  condition  that 
<v«,0>=0  (since  r  =0  at  x  =  0)  the  ion  velocity  at  each  cell 
boundary  can  be  calculated  by  starting  with  v*,i  and 
continuing  to  the  wall.  Using  equation  (50)  for  the  ion 
velocity  yields  the  following  expression  for  the  positive 
ion  flux: 


T+  :  (n.,i*i  +  n.,i)vi/2  (51) 

In  deriving  the  negative  ion  momentum  equation,  it  is 
assumed  that  the  negative  ion  velocity  will  be  dominated  by 
its  thermal  velocity.  Like  Von  Engel,  therefore  the 
inertial  term  (^<v>/9x)  can  be  dropped.  Then  for  a 
quasi-steady  state,  the  negative  ion  momentum  equation 
reduces  to: 

v_<v.>  =  -eE/m-  -  (  kT  - /n -m  -  )  3n  -  /  9x  (52) 

Equation  (52)  can  be  rewritten  as: 

r_  =  - D - 5n - / 5x  -  En -  (53) 
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In  their  present  form  equations  (54a)-(54c)  can  be 
solved  for  the  electron  densities  in  the  next  time  step 
(j  +  1)  (given  the  densities  and  electric  field  values  from 
the  j,k  time  step.  However  before  doing  this,  it  will  be 
beneficial  to  examine  the  method  employed  by  S.  Rockwood  in 
his  solution  of  the  Boltzmann  transport  equation  for 
electron  scattering  (Reference  11). 
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Consider  the  form  of  equations  (54a)-(54c)  if  the 
electric  field  and  n-  were  constants.  If  this  were  the 


case,  equations  (54a)-(54c)  could  be  written  in  the  form 
( 1  1 : 2349  )  : 


N 

n  i  =  Z  C i  , n  n  m  (55) 

M  =  1 

where  Ct,*  contains  the  coefficients  on  the  right  hand  side 
of  equations  (54a)-(54c)  (note  that  this  includes  the 
electric  field  and  n-;  since,  for  the  moment,  we  are 
assuming  they  are  constants).  If  we  now  assume  that  the 
densities  on  the  right  side  of  equation  (55)  are  in  the 
(j+1)  time  step,  it  is  possible  to  rewrite  equation  (55)  in 
the  form  of  a  matrix  equation  (11:2349): 

(I  -  T£)n(t+l)=  n ( t )  (56) 


where  1  is  the  identity  matrix.  One  advantage  of  equation 
(56)  is  we  can  now  solve  for  the  densities  implicitly.  A 
second  advantage  is  increased  communication  between  cells  is 
now  possible,  that  is  the  cells  are  more  effectively  coupled 
together.  These  N  coupled  equations  can  now  be  solved  by 
simple  matrix  methods. 

In  our  case  however,  the  electric  field  and  n-  are  not 
constants.  The  terms  involving  the  electric  field  and  n. 
are  non-linear  terms  (i.e.,  the  consist  of  products  of  two 
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densities  or  the  product  of  a  density  and  a  density- 
dependent  function).  Because  our  only  knowledge  of  the 
electric  field  and  n-  is  from  the  previous  tine  step,  the 
terms  involving  these  values  must  remain  in  the  j*fc  time 
step  (and  therefore  on  the  right  hand  side  of  equation  (56). 
The  modified  equation  (56)  is  (11:2350): 


(1  -  :C)  n  (  t  + 1 )  =  (1  ♦  ID  n  ( t  ) 


(57) 


where 

£’  is  a  constant  matrix  containing  the  coefficients  of 
the  linear  terms  of  equations  (54a)-(54c)  and  X  is  a 
tridiagonal  matrix  consisting  of  the  non-linear  terms. 
Since  all  of  the  values  on  the  right  hand  side  of  equation 
(57)  are  known  from  the  previous  time  step,  the  product  of 
the  right  hand  side  of  equation  (57)  is  a  vector  (which  we 
shall  denote  as  b).  Letting  &  -  ( 1  -  ' )  the  equation 

(57)  can  be  rewritten  as: 


& 


A  i  n  (  t  ♦  7  )  -  b 


(  58  ) 


Lets  look  now  at  equations  (54a)-(55c)  after  employing 
Rorkwood 1 s  solution  method: 


n  t  ,  \  (  1  —  X  (  z-3  )+DiI/hJ  )-(DiI/h*  )ni  .*  =  n't  ,  i  (  1 ♦d ini  ,  i  ) 


(  I  ^  /  2  h  )  E  >  (  n  tJ ,  j  ♦  n't  ,  j  ) 


i  =  1 


(5  9a 
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IW  I  I 


;  *  1  j  «  1  S*'. 

ni  (  1  - 1  (  z-o  ) ♦ 2D i I/h  *  )  -  <  D*  l/h  *  )  Ini  ,  .  *  i*n(  .  ,  .  i  )  = 

ni  ,  i  (  Udini  ,  i  )-<I*  /2h>  (E",  ( n  i  .  t  .  i+ni  .  .  )- 


E^i  - 1  ( n  i  ,  i  +  n  i  ,  j  - 1  )  ) 


2  s  i  s  N- 1 


j*i 


I  59b ) 


n  iV*  <  l-I(z-a)*DfI/h*+'tVT«/h)-<Df':/h,)n''c  ,  *  -  i  =  ni  ,»(  1  *d  Ln  -  ,  »  ) 

-(  „  i/2h)E*(nJ«  ,»  +  nJE  ,»-i  )  i=N  ( 59c  > 

Similarly  for  the  negative  ions: 

n'V }  i  (  1  -  to  )+D-T/h*  )-(D-X/h  *  Jni^js  ni  .  i  (  1  -dlni  .  i  ) 

-  ( ./2h  )E?  (ni.  .  »*ni  ,  * )  i  =  l  (60a) 


ni*,3t  (  1- to  )  ♦2D-l/h*  )  -D-t/h*  (n3/,1 ,  ♦  i+niVi  -  i  )  =  n-1  .  .  <  1  -din  i  ,  ,  ) 
ii  (  ni  ,  t  «  i  +  ni  ii)_  Ei-i(n-,i+n-,i-))] 

2  s  i  s  N-l  (60b) 


- ( Ic. /2h  )  [ E 


ni*,1!.  (  1  -  I?  )  ♦D-I/h  *♦  ( iv t  «/4h  )  )-  (D-I/h*  )n"-+.'  n  -  i=  ni . *  < 1 -din i  ,  »  ) 

- < _ _ I/2h )Ej < ni  ,  »*rii  ,  *  -  i  )  i  =  N  (60c) 

Substitution  of  the  positive  ion  flux  expression 
(equation  (51))  into  the  positive  ion  continuity  equation 
( equat ion  ( 8) )  yields : 


£ 
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ni,' »  =  n  ^  ,  ,  +  z  in  i  .  » - ( T/h )  Ini  ,  ,  v  i  -  n  ♦  ,  t  -  i  }  2  s  i  s  N-l 


(61a) 


Where  v ;  is  the  positive  ion  drift  velocity  given  by 
equation  (50).  Applying  the  boundary  conditions  to  obtain 


a 


4  1 


3 


n  «  , i  and  n  ♦  ,  *  we  get : 


n  ♦  *  i  =  n  ♦  ,  i  ♦ztni  ,  i  -  <  I/h  )n'l  ,  i  v{  1  =  1  (61b) 

n«,K=  n^,M*2ir/t.N-(I/h)(n*.iiVss  -  n  •  .  n  -  i  v  »  .  j  )  (61e) 

where  v*«  is  the  ion  sound  speed  (equation  (44)).  With  the 

exception  of  the  zmt  term,  equations  (61a)-(61c)  are  made 

s 

up  of  only  non-linear  terms  (because  is  a  function  of 
density).  Therefore  we  can  not  apply  Rockwood ’ s  method  for 
calculating  the  the  positive  ion  densities.  All  of  the 
values  on  the  right  hand  side  of  equations  (61a)-(6lc)  are 
known  from  the  previous  time  step.  Therefore  it  is  possible 
to  calculate  the  positive  ion  density  by  direct 
subs t i tut i on . 

Code  Operation 

Let’s  look  now  at  how  the  equations  developed  above  are 
used  in  a  computer  code  to  calculate  the  charged  particle 
distributions  and  their  associated  electric  field.  Looking 
at  the  basic  algorithm  will  make  it  easier  to  understand  the 
operation  of  the  code.  The  basic  steps  of  the  algorithm 
are  : 

1)  Input  the  initial  densities  and  electric  field 
(usually  based  on  Schottky  theory). 
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2)  Perform  a  LI'  decomposition  on  the  "A  matrix  of  both 
species  (i.e.  break  the  A  matrix  into  Upper  and  Lower 
Diagonal  matrices)  . 

3)  Calculate  the  "b"  vector  of  the  electrons  and  negative 
ions  and  solve  for  their  densities  by  using  back 

subs t i t  ut i on . 

4)  Calculate  the  positive  ion  velocity  for  each  cell 
boundary . 

5)  Using  the  velocities,  find  the  positive  ion  density 
f or  each  cel 1 . 

6)  Using  the  charged  particle  densities,  calculate  the 
electric  field. 

7)  Return  to  step  3.  Repeat  this  process  until  the 
solution  converges  to  steady  state. 

Let’s  now  look  at  each  of  these  steps  in  more  detail. 

As  we  see,  the  first  step  is  to  input  the  initial  densities 
and  electric  field.  These  are  the  values  the  code  will  use 
as  the  old  values  (i.e.  zeroth  time  step)  in  the  first  time 
step.  For  our  runs  in  the  two  component  case  the  amb l po lar 
profiles  were  used  as  the  initial  parameters. 

To  solve  the  matrix  equations  a  LU  decomposition  is 
performed  on  the  A  matrices  of  both  species.  Since  these 
matrices  are  constants,  the  LU  decompositions  must  be 
performed  only  once  by  the  code.  The  LU  decomposition  is 
performed  by  a  subroutine  called  LUDECOMP  taken  from 
Reference  10.  With  these  decompositions  completed  th*-  code 
begins  the  first  time  step. 
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Using  the  densities  and  electric  fields  from  the 
previous  time  step  (the  initial  values  for  the  first  time 
step),  the  b  vectors  are  calculated.  The  densities  are  then 
solved  for  using  back  substitution.  This  function  is 
performed  by  a  subroutine  called  Ll'BKSB,  also  taken  from 
Reference  10.  (This  is  actually  a  two  step  method  in  that 
the  operation  is  performed  first  on  the  electrons  and  then 
on  the  negative  ions). 

As  mentioned  in  the  development  section,  the  positive 
ion  equation  is  composed  only  of  non-linear  terms  so  it  is 
impossible  to  apply  the  matrix  method  to  solve  for  its 
densities.  To  calculate  the  positive  ion  densities,  the 
code  must  first  calculate  the  ion  velocity  for  each  cell 
boundary.  To  do  this,  the  code  uses  a  DO  loop  and  equation 
(50).  Because  the  velocity  at  the  center  (v0)  is  zero, 
equation  (50)  can  be  used  as  a  recursion  relation  for  the 
positive  ion  velocity.  The  velocities  for  all  cell 
boundaries  up  to  vN_,  (the  left  boundary  of  the  last  cell) 
are  calculated  in  this  way.  The  wall  velocitv  v  N  is  taken 
care  of  by  the  wall  flux  boundary  condition.  With  these 
velocities  calculated  the  code  then  uses  equations  (61a)- 
(61c)  to  calculate  the  cell  densities.  Like  the  velocity- 
calculations  this  is  done  on  a  cell-by-cell  basis  using  a  DO 
1  oop . 

Now  that  the  code  has  calculated  all  of  the  new  charged 
particle  densities,  the  updated  electric  field  can  be 
calculated.  As  in  the  illustrative  example  in  the  beginning 
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of  this  section,  the  electric  field  is  calculated  using  a 
expression  derived  from  Poisson’s  equation.  In  this  case, 
the  electric  field  expression  is: 

•J  ,  .  J  ,  1  j  J 

Ei=  ( e/ 1  o  )  i  ( n ♦  ,  m-  nt,M~  n - , h  )  h 
M  =  1 

The  calculation  of  the  electric  field  marks  the  end  of 
the  time  step.  At  this  point  the  code  increments  the  time 
and  returns  to  the  third  step  (calculation  of  the  b 
vectors).  This  procedure  is  repeated  until  the  density 
profiles  reach  steady  state.  The  following  criteria  will  be 
used  to  determine  when  the  profiles  have  reached  the  steady 
state.  We  shall  assume  that  a  steady  state  exists  when  the 
positive  ion  density  changes  less  than  IX  over  a  time  span 
equal  to  ten  characteristic  diffusion  times.  For  our 
preliminary  runs,  we  will  limit  the  time  iterations  by 
establishing  a  time  loop  counter.  When  the  loop  counter 
reaches  a  set  value  the  program  will  terminate.  The  purpose 
for  doing  this  it  to  allow  us  to  monitor  the  code  operation 
in  the  early  "checkout  phase"  of  the  code.  At  run 
termination  the  code  writes  the  charged  particle  densities 
and  electric  field  to  an  output  file  called  "fdjf.res".  A 
complete  listing  of  the  code  is  included  in  Appendix  C.  An 
outline  of  the  input  data  is  also  covered  in  Appendix  C. 


Results  and  Discussion 
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Up  to  this  point  we  have  only  been  concerned  with  the 
development  of  the  model.  We  want  to  now  turn  our  attention 
to  whether  or  not  this  model  correctly  predicts  the  behavior 
of  a  real  plasma.  Our  initial  plans  called  for  testing  the 
code  in  phases:  1)  one  component  diffusion;  2)  two 
component  plasma;  3)  three  component  (one  dimensional  slab); 
4)  three  component  (cylindrical).  The  purpose  for  dividing 
the  code  testing  into  four  phases  is  to  allow  the  code 
behavior  to  be  monitored  and  refined  in  problems  of 
increasing  complexity.  Unfortunately,  computational 
difficulties  arose  while  testing  the  code  in  the  two 
component  phase;  therefore,  the  discussion  in  this  chapter 
will  be  limited  to  the  first  two  phases.  The  third  and 
fourth  phases  will  be  discussed  in  the  final  chapter  of  this 
paper . 

Phase  One 

Phase  one  is  the  testing  of  the  code  in  a  one 
dimensional  diffusion  problem.  Our  system  in  this  problem 
consists  of  of  one  species  of  neutral  particles.  This 
simple  case  will  allow  us  to  see  if  the  code  is  capable  of 
modeling  a  diffusion  problem  in  which  only  density  driven 
diffusion  exists  (ie.  no  electric  field).  It  will  be 
assumed  that  the  particles  obey  the  continuity  equation  (37) 
and  flux  equation  equation  (38)  (minus  the  electric  field 


terms);  therefore,  the  code  will  be  calculating  solutions  of 
the  equation: 

5n/^t  =  Di(e?n/9t*)  +  zn  (62) 

To  determine  the  behavior  of  the  model,  this  one  component 
system  will  be  tested  for  three  different  cases.  These  case 
are:  1)  no  production  and  leakage;  2)  non-zero  leakage  and 
no  production;  3)  non-zero  leakage  and  production.  For  each 
of  these  cases,  we  shall  assume  that  the  initial  density 
profile  is  a  cosine  function  with  a  maximum  at  the  center 
(x=0)  and  a  minimum  at  the  wall.  To  determine  if  the  code 
is  properly  modeling  the  problem,  the  code’s  results  will  be 
compared  to  analytic  solutions. 

Lets  look  briefly  at  the  expected  results  of  the  three 
cases.  For  case  one  (no  production  or  leakage),  the  steady 
state  density  profile  is  apparent.  Because  no  production  or 
losses  occur  in  this  case,  the  density  profile  will  become 
uniform  in  time  (i.e.,  the  density  in  each  cell  will  be  the 
same ) . 

For  case  two  (losses  with  no  production),  lets  look  at 
the  solution  to  equation  (62)  derived  in  Appendix  B. 

n(x,t)=  noexpl  (z-D- ’it  ]cos(  >x)  (B-10) 

where  >  is  the  root  of  equation  ( B  —  7 )  (i.e.,  an  eigenvalue). 

The  loss  in  this  case  is  due  to  removal  of  particles  at  the 


boundary  wall. 


This  loss  is  given  by:  (n(L)v/4),  where  v  is 


the  velocity  of  the  particles  at  the  wall.  In  this  case  the 
production  term  (z)  is  zero;  therefore  equation  (B-10) 
predicts  that  the  density  profile  will  be  a  cosine  function 
decaying  in  time. 

We  again  turn  to  equation  (B-10)  in  evaluating  the 
third  case.  In  this  case  we  will  choose  a  z  that  makes  the 
exponent  in  equation  (B-10)  zero,  i.e.,  the  steady  state 
solution.  From  this  case  we  will  be  able  to  determine  if 
the  code  is  capable  of  reaching  a  steady  state  solution  in  a 
problem  in  which  particles  are  being  both  lost  and  created. 

The  code  produced  the  expected  results  for  each  of  the 
three  cases.  Since  the  profiles  of  cases  one  and  two  are 
uninteresting  (i.e.,  flat  and  zero,  respectively),  we  will 
only  discuss  the  case  three  results  in  detail.  A  comparison 
of  the  code  predicted  profile  and  the  analytic  solution  is 
shown  in  Figure  7.  Plotted  is  the  normalized  density  versus 
the  cell  number  for  a  steady  state  solution  (with 
De=  5.17x105  cm2/sec;  z-  1.195x10®  sec*1;  [ eq .  (B-7),  L=l] 

-=  1.517).  The  analytic  solution  densities  were  calculated 
at  the  right  hand  side  of  the  cells  (i.e.,  at 
x  =  . 1 , . 2 , . 3 , . . . , 1  ).  As  can  be  seen  from  Figure  7,  the 

code  and  solution  results  are  in  close  agreement.  The 
largest  difference  occurs  in  cell  7,  with  a  percent 


difference  of  4 . 3%  . 


Code  Predicted  Density  vs  finalytic  Solution  Dens 


An  unexpected  result  of  this  phase  was  that  of  code 
sensitivity  to  round-off  error.  Attempts  to  run  the  model 
with  single  precision  variables  resulted  in  crashes  of  the 
computer  code.  This  problem  was  corrected  by  performing  the 
calculations  with  double  precision  variables. 

Phase  Two 

Having  successfully  completed  the  first  phase,  testing 
of  the  code  for  the  two  component  plasma  began.  The  purpose 
of  this  step  is  to  increase  the  complexity  of  the  problem 
and  test  the  behavior  of  the  code.  As  phase  one,  we  have 
only  one  loss  and  source  process;  however,  we  now  have  the 
additional  features  of  an  electric  field-driven  flux  and 
coupling  of  two  oppositely  charged  species  by  the  electric 
field. 

One  item  of  particular  interest  in  this  phase  will  be 
the  effect  of  the  electric  field  on  numerical  instability. 
The  addition  of  the  electric  field  will  make  the  plasma 
equations  (finite  differenced  continuity  and  momentum 
equations)  sensitive  to  the  choice  of  time  step.  Bv  looking 
at  the  time  scales  of  the  electrons  and  ions  (or  more 
appropriately  the  difference  between  the  time  scales)  it  is 
possible  to  understand  how  the  electric  field  sensitizes  the 
plasma  equations.  The  particle’s  time  scale  is  the 
characteristic  time  in  which  the  species  reacts  to  changes 
in  the  electric  field  and  density  gradients.  The  electron 
time  scale  is  short  compared  to  that  of  the  ions.  Therefore 
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during  a  large  lime  step,  the  electron  density  can  change 
greatly  over  that  of  the  ions.  This,  in  turn,  can  cause 

y 

large  changes  in  the  electric  field.  These  field  changes 
are  then  fed  back  into  the  density  calculations  of  the  next 
time  step.  The  overall  results  of  a  large  time  step  can  be 
erratic  density  and  electric  field  profiles. 

As  in  the  first  phase  the  code  results  will  be  compared 
with  another  solution  to  determine  if  the  model/code  is 
correctly  modeling  the  plasma.  We  will  use  Self  and  Ewald’s 
model  as  the  comparison  solution  in  this  phase.  Self  and 
Ewald’s  model  is  ideal  for  this  comparison  for  three 
reasons.  The  first  reason  is  that  Self  and  Ewald’s 
solution,  like  our  model,  is  in  slab  geometry.  Therefore  a 
more  direct  comparison  can  be  made  using  this  model.  The 

i 

r  second  reason  is  that  the  boundary  conditions  applied  by 

Self  and  Ewald  are  similar  to  ours,  i.e.,  the  Bohm  criterion 
is  applied  for  the  ions  at  the  wall.  The  third  reason  is 
that  the  results  of  Self  and  Ewald’s  model  can  be  easily 
compared  for  different  ionization  rates.  This  makes  it 
possible  to  test  our  model  in  different  cases,  as  was  done 
in  the  first  phase. 

In  phase  one  we  simply  varied  the  source  rate  to 
develop  the  different  cases  in  which  the  code  was  tested. 

In  this  phase  we  will  indirectly  vary  the  source  rate 
(ionization)  by  testing  the  code  in  different  background  gas 
pressures  (the  ionization  rate  is  proportional  to  the 

%  neutral  density  which  is  dependent  on  the  gas  pressure). 
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This  can  best  be  understood  by  a  quick  review  of  Self  and 
Ewald’s  solution.  The  density  profile  predicted  by  Self  and 
Ewald’s  model  is  strongly  dependent  on  the  pressure  of  the 
background  gas.  This  can  be  easily  seen  in  Figure  2,  where 
the  density  profiles  for  several  different  pressures  are 
plotted.  In  this  plot  the  pressure  dependence  is 
represented  by  the  pressure  parameter  (A)  (14:2487).  A  is 

the  ratio  of  the  ionization  rate  to  the  sum  of  the  ioniztion 
rate  (z)  and  ion-neutral  collision  frequency  ( 0  ).  For  low 
pressure  gases  the  ion-neutral  collision  frequency  is  small 
so  the  pressure  parameter  is  approximately  unity.  Therefore 
a  low  pressure  gas  is  a  gas  in  which  the  mean  free  path  for 
collisions  between  ions  and  neutrals  is  large  compared  to 
the  discharge  tube  radius  (14:2487).  Just  the  opposite  is 
true  in  a  high  pressure  gas.  In  a  high  pressure  gas  A 
approaches  zero  and  the  mean  free  path  of  ion-neutral 
collisions  is  much  greater  than  the  tube  radius. 

The  three  particular  cases  we  will  look  at  are  A--1 
(  ..  <  <  1  )  ,  A=  .  5  (  v  -z  )  ,  A'tO  (  ,  >  >  z  )  .  In  the  low  pressure  case 
the  density  profile  should  be  approximately  that  of  the 
profile  predicted  by  the  free  fall  theory  (14:2486),  see 
Figure  2.  The  second  case  is  an  intermediate  pressure  and 
the  expected  profile  lies  between  that  of  the  free  fall  and 
ambipolar  profiles.  The  last  case,  A-0,  is  a  high  pressure 
case.  The  expected  profile  is  like  that  of  the  ambipolar 


As  was  expected,  the  code  was  sensitive  to  the 


magnitude  of  the  time  step  interval.  It  was  found  that  when 
using  ten  cells  to  model  the  plasma,  the  time  step  had  to 
less  than  or  equal  to  10**  seconds  for  the  code  to  remain 
numerically  stable.  When  using  100  cells  the  time  interval 
was  restricted  to  less  than  10- >°.  These  restraints  forced 
long  computational  times  (e.g.,  to  bring  the  real  time 
calculation  of  the  ten  cell  model  to  one  millisecond 
required  100,000  time  steps). 

In  all  attempts  to  calculate  the  densities  of  the  two 
component  plasma,  the  positive  ion  density  became  ionization 
as  the  integration  in  time  progressed.  This  is  best 
demonstrated  by  the  positive  ion  profile  evolution  depicted 
in  Figure  8.  The  resulting  electric  field  profile  is  shown 
in  Figure  9. 

From  Figure  8  we  see  that  the  ion  density  becomes 
increasingly  erratic  as  the  integration  in  time  progresses. 
We  note  that  an  oscillation  first  becomes  apparent  at 
50„sec.  This  oscillation  then  spreads  inward  toward  the 
axis.  Allowing  the  calculations  to  progress  further  in  time 
resulted  in  negative  densities  being  predicted;  therefore, 
we  are  assured  that  this  behavior  is  code  generated  and  not 
rea  1  . 

The  cause  of  the  numerical  instability  (oscillations) 
was  not  discovered  during  the  course  of  the  research. 

During  the  writing  of  this  report,  code  errors  were 
discovered  in  the  ion  velocity  calculations.  These  errors 
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Positive  Ion  Dens 


were  at  least  in  part  the  cause  for  the  oscillations  se^n  in 
the  code  results.  A  detailed  analysis  of  the  results  of  the 
corrected  code  was  not  possible  due  to  time  limitations.  A 
quick  examination  of  the  ion  velocities  revealed  that  a 
inconsistency  existed  at  the  discharge  wall.  The  velocity 
on  the  left  edge  of  the  last  cell  (cell  boundary  before 
wall)  was  greater  than  the  ion  sound  speed  (by  a  factor  of 
four  in  one  case).  This  inconsistency  led  to  a  "pile-up"  in 
the  last  cell  of  the  grid,  which  in  turn  led  to  oscillations 
in  the  ion  densities.  This  behavior  implies  that  our 
boundary  condition  for  the  positive  ions  (i.e.,  the  ion 
velocity  at  the  wall  is  exactly  equal  to  the  ion  sound 


speed)  is  invalid. 


VI.  Conclusions  and  Recommendations 


In  this  chapter  we  shall  discuss  some  of  the  areas 
which  should  be  explored  in  an  attempt  to  correct  the 
current  problems  experienced  in  the  model  and  code.  We  will 
also  discuss  the  final  two  phases  of  the  code  testing,  which 
could  be  undertaken  once  the  current  problems  are  corrected 
and  phase  two  testing  completed. 

Model/Code  Corrections 

The  first  area  which  must  be  researched,  before  code 
testing  can  continue,  is  the  boundary  condition  on  the 
positive  flux  at  the  discharge  wall.  In  particular,  the 
effect  of  allowing  the  ion’s  wall  velocity  to  exceed  the  ion 
sound  speed  should  be  examined.  As  an  alternative  to  fixing 
the  positive  ion  velocity  at  the  wall,  to  the  ion  sound 
speed,  is  to  use  equation  (50)  to  calculate  the  ion  speed  at 
the  wall. 

The  second  area  which  should  be  examined  is  that  of  the 
sensitivity  of  the  code  to  the  time  step.  Because  the 
electric  field  coupling  is  the  primary  cause  of  this 
sensitivity,  it  may  be  possible  to  reduce  this  sensitivity 
by  by  introducing  a  perturbation  term  to  the  electric  field. 
The  new  form  of  the  electric  field  would  be: 


E(x,t)=  E  o ( x )  +  E  » ( x , t  ) 


where  Eo  is  a  time  constant  field  term,  based  on  the 
expected  electric  field  profile.  For  example  in  the  case  of 
the  two  component  plasma,  Eo  would  be  a  Tangent  function,  as 
predicted  by  the  ambipolar  theory.  The  perturbation  Ej 
would  be  the  difference  between  the  field  calculated  using 
Poisson’s  Equation  and  Eo.  Therefore  all  changes  in  the 
electric  field  from  one  time  step  to  the  other  would  be 
contained  in  the  perturbation  term.  We  would  expect  this 
perturbation  term  to  be  small;  therefore,  a  larger  time  step 
should  be  possible. 

This  change  will  also  better  couple  the  electric  field 
into  the  density  calculations.  This  coupling  can  be  seen  by- 
substituting  this  new  electric  field  expression  in  the 
electron  and  negative  ion  density  equations  (equations 
((59a)-(59c)  and  ( 60a ) - ( 60c ) ) .  Currently  the  electric  field 
terms  are  retained  on  the  right  hand  side  of  these  equations 
(because  the  are  non-linear  terms).  By  writing  the  electric 
field  with  a  constant  and  a  perturbation  it  will  be  possible 
to  bring  the  E0  terms  to  the  left  hand  side  of  the 
equations.  Substituting  for  the  electric  field  yields  the 
following  for  the  electrons: 

niVi  ( 1 -T( z-o )♦ ( DtT/h* )  +  ( * ft/2h )E. . i )- { ( Dii/h  * )♦ 

(pjl/2h)E8,ilni,i=  n  £  ,  i  (  1  ♦d  in'-  ,  i  )  -  (  *,  ^  t  /  2h  )  {  Ei  ,  ■ 

(  n  i  ,  i  +  n  i  ,  j  ) )  i  r  1 


(  63a  ) 


nt  , ,  {l-l(z-o)  +  (2Dil/h* )*(«ti/2h)(E,,t-E0l|.il|. 


ni  ,  i  -  i(  (  x.^  I/2h  ]  E  o  ,  t  -  i  +  (Dti/hl]  )  ♦  n  e  ,  »  ♦  1  (  [  l  /  2  h  ]  E  o  ,  ■  ♦  i  ■ 


IDiI/h*]  >=nl  .  i  (  l+dxnJ.  .  ,  )-(  u  I/2h)(EJ,  .  ,ni  ,  ,-E5,  .  ,ni  .  i  -  ,  ) 


i  <  i  s  N  - 1 


(63b) 


n{,«  1  1  -  X  (  z-o  )  +  (DtI/h*  )  +  (  u  l/2h  >E0  .  n+Tvt.  .  t/h)  +nJt  .  N  -  i  <  I  -  I/2h  ] 


l  Di  I/h  2  ]  )  =n  t  ,  n  (  1+dlnl  ,  h  )  -  (  u  EI/2h  )  E  }  ,  N  (  n  2  ,  h+t.  i  ,  N  -  i  ) 


i  =  N 


(  63c  ) 


Similarly  for  the  negative  ions: 


n-,i[l-Ia+(Dil/h»)  +  U  I/2h)E«  ,  i  )-(  ( D-T/h  2  )  ♦ 


{.»  I / 2 h  ) E o  ,  i  } nJ-+  :  n -  .  » (  1  -dln£  .  i  )  -  (  „  I/2hM  E  ,J ,  , 


(  n  ,  j  +  nr,,)) 


(64a) 


ni  V,  {  l-Ic+(  2D-i/h  *)  +  (..  T/2h  )  (E»  ,  ,-Eo 


n'-Vi  -  i  (  [  -  .  l/2h]E,  ,  .  -  ,  +  (D-l/hM  )+n -*;  ,  ,  !  (  U.  '/2h  )Eo  .  ,  ♦  . 


(D-~/h2]  )=n-l  ,  i  (  1  -d  in  i  .  (  )-(^.l/2h)(E'i  ,  in  -  ,  ,  .  -  E  -  i  n".  ,  t  -  ) 


is  i  s  N  -  1 


(  64b  ) 


nJ-  ,  n  (  (  D-  I/h  *  I  +  (  „  :/2h)Eo,«  +  lvTi,./h|+n'.,\.,(U  l/2h  1- 


[  D-i/h  *  )  )  =n:l  ,  h  (  1-dini  .  »  )-  (  t  .l/2h  )E  j  ,  »  (  n :  ,  N  +  nJ-  ,  N  -  ,  ) 


i  =  N 


(  64c  ) 


*  j  i  1  j  i 

where  n‘i(  m:  (ni,i+  n  *  ,  »  ♦  j  )  ;  n«,i-=  (  n  •  .  ,  ♦  n*f»-i) 


Once  the  two  areas  discussed  above  have  been  addressed, 
it  will  possible  to  continue  the  testing  of  the  code,  as 
outlined  in  phase  two.  In  this  section  we  will  discuss  the 
code  testing  which  should  occur  after  the  phase  two  testing 
is  complete. 

Having  successfully  completed  phase  two,  we  will  have 
proved  the  finite  differencing  method’s  ability  to  model  a 
simple  plasma.  Therefore  our  primary  goal  in  the  third 
phase  is  to  establish  the  general  density  profiles  of  the 
three  component  plasma  for  different  sets  of  parameters 
(i.e.,  different  ionization,  attachment,  and  detachment 
rates).  A  secondary  goal  will  be  to  compare  the  code’s 
results  with  Lee's  and  Von  Engel’s  solutions  to  verify  the 
relationships  among  the  particle  densities  (i.e.,  greater 
positive  ion  densities  than  electron  densities  throughout 
the  volume,  negative  ions  confined  to  the  axis  region,  etc.) 
reported  in  their  papers. 

Though  we  cannot  expect  profiles  identical  to  Von 
Engel’s  or  Lee’s  since  their  results  are  in  cylindrical 
geometry,  a  comparision  of  our  model’s  results  with  their 
results  should  be  possible.  The  basic  relationships  between 
the  profiles  of  the  different  species  should  be  comparison 
in  both  geometries.  From  this  comparison  it  may  be  possible 
to  draw  conclusions  on  the  accuracy  or  inaccuracy  of  their 
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APPENDIX  A 


Development  of  Self  and  Ewald's  Plasma  Model 

The  authors  start  with  the  following  steady  state 
continuity  and  momentum  equations: 


t  (nv|:G  ( A-  1  ) 

n(v  v  )  v  +  Gv:  (  nq/m  )  E- (  7  *r/m)-n.v  (A- 2) 

where : 

v  is  the  collision  frequency  for  momentum  transfer; 

G  is  the  volumetric  production  rate; 

>r  is  the  pressure  tensor. 

Making  the  following  substitutions: 

v  tr  =  kTrn  (  A-  3  ) 

G  =  o  j  n  £  (  A  -  4  ) 

Where  -j  is  the  ionization  frequency.  Equation  ( A  —  1  )  can 
be  rewritten  by  substituting  (A-4)  into  ( A  —  1  )  : 


v  (  7  n  )  +  n  (  7  v  )  =  ^  ,  n 

v  (  v  n  )  (  1  /n  )  +  v  v  -  u  , 

v  c l Ln ( n ) ]  +  ( r  v ) =  0  j  (A-5) 

Substituting  <  A  —  3 )  in  ( A  —  2 )  yields: 


n(v  7 ) v  +  n.  v  =  nq ( E/m ) - ( kT/m ) 7n  -  n.  v 
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dividing  by  n  and  letting  o'  =  oj  ♦  j  gives 


(v  y)v  =  q(E/m)-(kT/m)y{Ln(n ) )-  o'v  (A-6) 

In  one  dimensional  slab  geometry  (discharge  between  to 
infinite  parallel  plates,  equation  (A-5)  reduces  to: 

dv/dr  +  v  d(Ln(n)]/dr  =  v  j  (A-7) 

where  r  is  the  spatial  variable.  Substituting  E=  -ye  for 
the  electric  field  and  -e  for  q  in  equation  (A-6)  yields  the 
momentum  transfer  equation  for  the  electrons: 

v  dv/dr  =(e/m£)dc/dr  -(kT£/m£)  d(Ln(n)]/dr  -  u£v  ( A- 8  ) 

0 

similarly  for  the  ions 

v  dv/dr  =- ( e/m  ♦  )d: /dr  -(kT./m,)  d[Ln(n)]/dr  -  .)  v  (A-9) 

Note  that  the  author  assume  quasi-neutralitv  and  therefore 
make  no  distinction  in  the  density  variable  (i.e.,  n=n£=n«). 

Equations  (A-7)-  (A-9)  are  normalized  using  the 
following  factors. 

r=e«/kTi  ;  u=  ( m  ♦  /kT  £  )  1  ■'  *  v  ;  s=(m,/kTi)1  *  oj  r  ; 

R  ’  =  (  m  ,  /m  t  )(.,/'  e  )  *  A=  o  ]  /  v+  :  I  =  T  ♦  /  T  £  (  A  —  1  0  ) 
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Using  the  chain  rule  yields  an  expression  for  d  /dr 


d/dr  =  (  m  . /kT  i  )  1  '  *  d/ds  (  A  -  1  1  ) 

Direct  substitution  of  ( A  —  1 1 )  and  ( A -  1 0 )  in  equation  ( A  —  7  ) 
yields  the  normalized  continuity  equation: 

d(u)/ds  +u  d[Ln(n)]/ds  r  A  ( A- 1 2  ) 


Before  normalizing  the  electron  momentum  transfer  equations 
the  authors  assume  that  the  electron  inertial  term  can  be 
ignored  due  to  the  small  mass  of  the  electron  (14:2487). 

Now  performing  the  normalization  gives: 


u=R’(d'/ds)-  I  dlLn(n)]/ds 


(  A-  1  3  ) 


Normalization  of  the  ion  equations  yields: 

u+u(du/ds  I  :  -d'./ds  -  "  d[Ln(n)]/ds  (  A  -  1  4  ) 

Th«  authors'  next  step  is  to  make  u  t  hf  independent 
variable  in  equat  ions  ( A- 1 2 ) - ( A- 1  4  )  .  The  chain  rule  is  used 
to  obtain  an  expression  for  d/ds.  This  expression  is 
substituted  into  equations  (A-12)-(A-14)  yielding: 

1  ♦  u  d[Ln(n)l/ds  =  A  ds/du  ( A -  1 5  ) 

u  ds/du  =  R’(dr/du  -  d|  l.n(n)  ]/du)  (A -16) 
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u( 1  +  ds/su  )  =  -d"/du 


I  d ( Ln ( n ) ] /du 


(  A-l 7 ) 


Equations  <  A  —  1 5 ) —  ( A  —  1 7 )  are  now  solved  to  find  ds/su, 
s,  Ln(n/n«),  and  r.  (  no  is  the  density  at  x  =  0).  Solving 
(A- 16)  for  dri/ds  we  find: 


dr /ds  =  (u/R')ds/du  +  d(Ln(n))/du  ( A- 1 8 ) 

Substituting  ( A- 1 8 )  for  dr/ds  in  (A-l?)  yields: 

d(Ln(n)]/du  =  (A/u)ds/du  ( A  —  1 9 ) 

(A-19)  can  now  be  used  in  (A-15)  to  solve  for  ds/du  . 

ds/du  =  [  (  l  +  D-u*  ]  /  l  A  (  l  +  I)  +  <  1+1/R’  )u*  )  (A-20) 

( A  —  2 0 )  is  integrated  to  find  s. 

s=[(l  +  *.  )‘  *1  Ul/R’+A)  ]  /  l  A  »  M  1  +  1/R')  3  3  ]Tan*  MWu)-u/(  1  +  1/R'  ) 

(  A  -  2  1  > 


where  V=l  (  1  +  1/  R  ’  )  / <  A (  l  +  I)  )  *  3  1  . 

Integration  of  (A-19)  (note:  n ( 0 ) =  0  )  yields 


Ln(n/n» )=-[ 1+1/R '+Al/(2( 1+1/R ') lLn(Y) 


(  A-22  ) 
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APPENDIX  B 

Development  of  the  Solution  for  One  Component  Diffusion 

This  appendix  outlines  the  development  of  the  solution 
for  one  species  diffusion  in  one  dimension.  The  system 
consists  of  two  infinite  parallel  plates  separated  a 
distance  2L  .  It  is  assumed  that  the  particles  conservation 
is  governed  by  the  following  continuity  equation: 


5n/ St  =  - ( c  T)  +  zn 


(  B  —  2  > 


where  z  is  the  particle  production  rate.  It  is  further 
assumed  that  the  particle  flux  can  be  expressed  by  Fick’s 
Law : 


r  =  -Dr  n 


(  B-2  ) 


where  D  is  a  constant  diffusion  coefficient. 

As  a  boundary  condition,  it  is  assumed  that  the  flux  at 
the  wall  (i.e.,  the  particle  loss)  is  given  by: 

r  (  L  )  =  (  n  (  L  )  v  )  /  4 

where  v  is  the  average  velocity  of  the  particles. 
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Casting  equations  <  B  —  1 )  and  <  B  —  2 )  in  cartesian 


coordinates  and  substituting  (B-2)  into  (B-l)  yields: 

3n/?t  =  D~?n/ox*  +  zn  (B-3) 

Applying  the  method  of  separation  of  variables  we  assume 
that  n(x,t)=X(x)T(t).  Substituting  for  n  in  equation  (B-3) 
gi  ves  the  following  expression: 

X ( x ) 3T ( t ) / 3t  =  DT(t)3?X(x)/3x*  ♦  X(x)T(t)z 


or 


(  1/D)  [TMtl/TIt)  1  -  (z/D)  =  X  "  (  x  )  /X  (  x  ) 

The  only  wav  for  this  equality  to  be  true  is  for  both  sides 
of  the  equation  to  be  equal  to  a  constant  (  ~>  )  .  Therefore  we 
get  to  separate  equations. 

X " ( x ) /X ( x  )  =  -  ( B- 4  ) 

(  1/D)  { (T ’ (t )/T(t)  ]  -  z)  =  >  (B-5 ) 

Expanding  (B-4)  we  get  the  following: 

X  *  (  x  )  =  >  X  (  x  ) 
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Before  we  can  solve  this  differential  equation  we  must 
establish  the  boundary  conditions  for  the  separated 
equations.  The  first  boundary  condition  comes  from  the 


boundary  condition  of  the  flux  at  x=0,  i.e.,  r(0,t)=0.  If 
the  flux  at  x=0  is  zero  then  -DT ( t ) X ’ ( x ) =0 .  Since  this  must 
hold  true  for  all  times  X'(0)=0.  Similarly  the  boundary 
condition  on  the  flux  at  x  =  R  yields  that  X  ’  ( R )  =  - ( v/4D ) X ( R ) . 

Now  we  can  begin  to  solve  the  differential  equations 
above.  The  first  step  is  to  determine  if  the  constant  1  is 
negative,  positive,  or  zero.  Assuming  first  that  1  is 
positive  and  that  >=>  .  Using  this  in  the  differential 

equation  we  find  that  the  solution  is: 

X(x)=  Aiexp(->.x)  ♦  A»exp(/x) 

Taking  the  derivative  of  this  solution  gives 

X‘(x)  =  -Ai-exp(-*x)  ♦  A  i  •  e  x  p  (  x  ) 

From  the  boundary  condition  at  x=0  we  find  that  Ai=Aj. 
Substituting  for  At  and  applying  the  boundary  condition  at 
x  =  1,  yields: 

Ai/[exp(->L)  -  e  x  p  (  *  L  )  ]  =  -<v/4D)Ai(exp(-.*L)  +  exp(  L)] 

Therefore,  - =- ( v/4D ) ( 1 /Tanh ( a R ) ) ;  however,  there  are  no 
roots  to  this  equation.  Therefore,  >  >0  is  not  a  possible 


solution. 


A  second  possible  solution  is  that  -=0.  Applying  this 
we  get : 


X  (  x  )  =  Ai\  ♦  A] 

and 

X  '  (  x  )  =  A, 


Applying  the  boundary  condition  at  the  x=0  tell  us  Aj  must 
be  zero.  If  Ai  =  0  then  the  boundary  condition  at  the  wall 
can  not  be  met . 

The  only  choice  left  is  that  -  is  negative.  Assuming 
yz-\  we  get: 


X ( x  )  =  AiCos(/x)  +  Aisinl  <\ ) 


(  B-6  ) 


X  ’  (  x  )  =  -Ai'Sin(>x)  +  A  *  •.  c  o  s  (  •  x  ) 


Because  X’(0)=0,  Aj  must  be  zero.  Using  this  and  applying 
the  wall  boundary  condition  yields: 


- A i  s i n ( • L )  =  -(v/4D)AiCos(-L> 

Therefore  *  is  given  by  v= ( v/4D ) cot (  • R  )  .  This  relation  has 
possible  roots  and  thus  is  our  solution. 

X(x)=  AicosUx)  where  • = < v/4D )cot ( • R )  (B-7) 
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Turning  our  attention  to  the  the  time  dependent 


The  solution  to  this  equation  is 

T ( t )  =  exp  <  [  z-D>?  ]t )  (B-8  ) 


Since  n ( x , t ) =X ( x ) T ( t )  we  get  for  our  finial  solution 


n(x,t)=  A i |exp( [ z-Da ? ) t )cos ( Xx  )  }  (B-9) 


equation  we  find: 

T  ’  ( t )  -  (z-D>?  )  T  ( t  )  =  0 


Defining  the  density  at  t=0  and  x=0  as  no  we  can  evaluate 
A i ,  i.e.,  n(0,0)=no»  so 


n(x,t)=n0{exp< [ z-D>  ?  J  t ) cos ( *  x  )  ) 


(  B-  1 0  ) 


where  -  is  the  roots  of  - •  =  ( v/4D ) cot ( * L ) .  It  should  be  noted 
that  this  solution  also  can  be  used  as  a  solution  to  the 
ambipolar  model  if  we  substitute  D*  for  D. 


APPENDIX  C 


Finite  Difference  Model  Code 

This  appendix  is  a  listing  of  the  latest  computer  code 
used  in  the  research.  The  program  requires  input 
from  a  data  file  called  fdif.dat.  The  data  in  fdif.dat  is 
as  f ol lows : 

ts  -  #  of  time  steps  to  perform 
c  -  #  of  cells  in  grid 

D#-  Electron  Diffusion  Coef. 
u  -  Electron  Mobility 

T*,T,-  The  electron  and  ion  temperatures  in  Kelvin 
.?  -  The  ion-neutral  collision  cross  section 
z  -  Ionization  rate 
m  ,  -  I  on  mass 

s  -  k/h  ratio  of  time  step  interval  to  cell  width 
?  -  Attachment  rate 
d  -  Detachment  rate 

ns-  *  of  species  in  plasma  (2  or  3) 

The  program  outputs  the  particle  densities  and  electric 
field  for  the  last  time  step  into  a  data  file  fdif.res. 


*  PROGRAM  FDIF.FOR  --  FORTRAN? 7  --  JOHN  LITAS 

*  THIS  PROGRAM  CALCl ALTES  THE  PARTICLE  DENSITIES  IN  A  TWO  OR 

*  THREE  COMPONENT  PLASMA. 

implicit  real*8  (a,e,o-z) 

real *8  nel( 10) ,  ne2  ( 10) ,npl ( 10) , np2 ( 1 0 ) 

real *8  el ( 0 : 1 0 ) ,nnl ( 10 ) , nn2 ( 1 0 ) , rho ( 10 ) 

real *8  v(0:10), kb, k,fe (0:10), fp (0:10) 

real*8  nO  , nu , 1 amda 

integer  ts 

real*8  ae ( 1 0 , 1 0 ) , an i ( 1 0 , 1 0 ) 
integer  indxn ( 1 0 ) , i ndxe ( 1 0 ) 

open (unit=10, files' fdif.dat’ ,status='old' ) 
open (unit=ll, files’ fdif. res’  ,status=’old’  ) 

read ( 1 0 , * ) ts 
read ( 10 , *  )c 
read ( 1 0 , * )de 
read ( 10 , * )be 
read ( 1 0 , * ) te , tp 
read ( 1 0 , t ) press 
read( 10, * ) sigma 
read ( 10 , *  )  z 
read ( 1 0 , * ) s 
read ( 1 0 , t ) pim 
read! 10, * )alpha 
read! 10, * )d 
read( 10, * )nsp 

*  Define  contstants  to  be  used 

h=  1  .  /c 
k  =  s*h 

kb= 1 . 3807E-23 
em  =  9 . 1095E-31 
ec= 1 . 6022E- 1 9 
pi  =  4  .  *atan ( 1 .  ) 
cl=100. *ec/8.8542E-12 
m= int ( c ) 

n0=press* 133.2*6.022E+23/(8.3144*tp*1.0E+06) 

cel=100.*sqrt(kb*te/(2.*pi*em) ) 

cpl=100.»sqrt(kb*tp/(2.>pi*pim)) 

cp2=100.*sqrt(kb*te/(2.*pi*pim)) 

lamda= l./(n0*sigma) 

nu=cpl/lamda 

n=m 

nep=m 

nnp=m 

*  Call  setup  to  get  initial  conditions 

call  setup ( e 1 , ne 1 , npl , nn 1 , nsp ) 


*  Set  up  Constant  Matrix  (Left  side)  for  electrons  and 

*  Negative  ion  (if  present) 

a j 1 =k  *de/h*  * 2 

aj2= ( 1 . -k* ( z-alpha )  +  a  j  1 ) 

a  j  3  =  (  1  .  -k  *  ( z-alpha )  +  2  .  *a j 1 ) 

ae ( 1 , 1 )  =a  j  2 
ae ( 1 , 2 ) = -a j 1 
do  5 , j j  =  2 , m- 1 
ae ( j j , j j- 1 ) =-aj 1 
ae ( j j , j j ) =a j  3 
ae ( j j , j j  +  1  )  =  -a j 1 
00005  continue 

ae(m,m-l )  =  —  a  j 1 
ae(m,m)=aj2+(k*cel/h) 

*  With  the  constant  matrix  calculated.  Call  the  sub.  Ll: 

*  Decomposition  to  find  the  decomp,  matrix  of  ae.  Decomp 

*  Matrix  will  become  the  new  ae  matrix. 

call  ludcmp(ae ,n,nep, indxe ,dae ) 

*  Repeat  procedure  if  negative  ions  are  present 

if  (nsp  .eq.  3)  then 
ajl=((k/h)**2./(l.+k*nu))*l kb* tp/pim ) 
a j2= ( 1 . +a j 1 ) 
a j3=  < 1  .  +2 . *a j 1 ) 

an i ( 1 , 1  ) =a j  2 
ani ( 1 , 2  )  =  -a j 1 

do  7 , j  j  =  1 , m- 1 
ani  ( j j  ,  j j-  1  )  =  -a j 1 
ani ( j j , j j ) =a j  3 
ani ( j j , j j+ 1 ) =-a j 1 
00007  continue 

ani ( m , m- 1 ) = -a j 1 
ani (m,m)=aj2+(k/h)*cpl 
call  1 udcrap ( an l , n , nnp , i ndxn , dn i  ) 
end  if 


*  Start  of  time  loop 

do  1 0 , j  = 1 , ts 

t2=be*k/(2.*h) 
t  3  =  ne 1 ( 1  )  +nel ( 2 ) 

1 1 =ne 1 ( 1 )*( 1  .  +k*d*nn 1 ( 1  )  ) 
ne2 (1  ) =t 1 +t2*t3*el ( 1  ) 


do  1 5 , i =2 , m- 1 
tl=nel ( i )*( 1 .  +k  *d*nn 1 <  i ) ) 
t3  =  e 1 ( i ) *  <  ne 1 ( i*l  )*nel ( i  )  ) 
t4=el ( i -1 )* (nel < i )+nel ( i-1 ) ) 
ne2 { i ) =t l  +  t2* ( t3-t4  ) 

00015  continue 

tl=nel(m)*(l.+k*d*nnl(m)  ) 
t3  =  el (m-1  1  * ( nel (m)+nel (m-1 ) ) 
ne2 I m )  =  1 1 - 1  2  * t  3 

call  i ubksb ( ae , n , nep , i ndxe , ne2 > 

00093  cont i nue 

*  Calculate  the  ion  drift  velocities 
00107  continue 
v  (  0  )  —  0  • 

vlp=2 . *h* lamda/ (2.*h+lamda) 
vt2=(nel (2)+nel ( 1 ))/(npl <2)*npl ( 1 ) ) 
vt3=vlp* vt2*z 
vt4=1.0E+04*ec*el(l)/pim 

v ( 1 )  =  . 5  * ( -vt 3  +  sqrt ( vt3**2 . +4 . *vlp*vt 4 )  ) 
ptl=npl(l)+z*k*ne2(l) 
pt2=(k/h)*v(l )*<npl C 1 )+npl (2) }/2. 
np2 ( 1 ) =pt 1 -pt2 

do  20 , i  =  2 , m- 1 

vt  2  =  (  ne  1  (  i  ♦  1  )+nel(i  )  )/<np2 ( i*l  )  -*-  n  p.  1  <  i  )  ) 
v  t  3  =  v  1  p  t  v  t  2  *  z 

v  1 4  =  1 .  0E  =  04*ec*el ( i  ) /pi m  +  v ( i-1  ) *  *  2 . / ( 2 . *h  ) 
v(i)=.5*(-vt3+sqrt(vt3**2.+4.*vlp*\t4)) 
pt]=v(i)*(npl(i+l)+npl(i)) 
pt2  =  v(i-l ) * ( n p 1 ( i-1 ) ♦np 1 (  i  )  ) 

np2(i  ) =np 1 ( i )  +  z*k*ne2<  l  )-(k/(2.*h)*(pt  1  - p t  2  )  ) 
00020  continue 

pt  1  :  .  5  *  v ( i-1  )  *  (  n  p  1 ( i-1  >  +  n  p  1 <  m  ) ) 
pt2=(k/h)» (npl (m)*cp2-pt  1  ) 
np2 ( m ) =npl (m)+z*k*ne2(m)-pt2 
00108  continue 

*  If  negative  ions  are  present  setup  b  matrix 

i f  ( nsp  . eq .  3 )  then 

tl=ec*k**2./(2.*h*pim*(l.+k*nu>) 
t  2  =  < 1 . -k*d*nel ( 1 ) ) *nnl ( 1 > 
t  3  =  e 1 ( 1 )* (nnl ( 1 )+nnl <  2  )  ) 
nn2( 1  )  =  t2  +  k*alpha*ne2(  1  )  ♦  1 1  *  t  3 
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do  25 , 1 =  2 , m-  1 
1 4  :  t  3 

t2=(  1 . -k*d*ne  1  (  i  1  )  *nn  1  (  i  ) 
t3  =  el(i)*(nnl(i)+nnl(i-M)) 
nn2(  i  )=t2+k*alpha*ne2(  i  )  ♦  t_  1  *  (  t3-t4  ) 
00025  continue 

t  2  = ( 1 . -k  *d  *  ne  1 (  m  )  1  *  nn  1 (  m  ) 
nn2(m)=t2+k*alpha*ne2!ro)-tl*t3 

call  1 ubk sb ( an i , n , nnp , i ndxn , nn 2  ) 
end  i  f 
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*  Calculate  the  Electric  field 


el  ( 0  )  =0 . 
do  30 ,  i  =  1  ,  m 

rhol i )=np2( i  )-ne2( l  )  - nn  2 (  i  ) 
el(i)=el(i-l  )+cl*h*rho(i  ) 

00030  continue 

*  Set  old  densities  equal  to  newly  calculated  densities 

*  for  next  time  step 

do  35  ,  i  =  1  ,  m 
ne 1 ( i  )  cne 2 (  i  ) 
npl ( i ) =  np2 I i ) 
nnl(i)=nn2(i) 

00035  continue 


00010  continue 
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t i me: k  * ( j  -  1  ) 

write(llt*)’h=',h,’k=’,k,’ 
wntelll,*)'  ’ 
write!  11, * )  ’ F 1 ec t  ron  Den s i 
write!  11  ,200)  <ne2<  i  ),ne2(  l 

♦  4  )  , 

X  1  :  1  ,  6 , 5  ) 

write! 1  1  ,*  )  ’  ’ 

write!  1  1  ,  *  )  'Positive  Ion  Density’ 

write!  1  1  , 200 ) ( np2 (  l  )  , np2 (  i ♦ 1  )  , np2 ( i ♦ 2 )  , np2 (  i  +  3 )  , np2 (  i  | 

+  4  )  ,  ? 

X  l  =  1  ,6,5)  \ 


t  i  me  :  ’  ,  1  i  me 

tv-  j 

♦  1  )  ,  n  e  2  (  i  ♦  2  )  ,  n  e  2  (  i  +  3  )  ,  n  e  2  (  i 


i  f  (  nsp  .  eq  .  3  )  t  hen 

w  r  i  t  e  (  1  1  ,  *  )  ’  ’  w] 

write!  1  1  ,*  I'  Sr-gnt  l  ve  I  on  Dens  1 1  y  ’  | 

writp(ll,200)(nn2(i),nn2(i  +  l),nn2(i-»2),nr;2(i+3),nn2( 


i  ♦  4  )  , 


X  i  =  1  ,  6 , 5  ) 
end  i  f 

write!  11,*)’  ’ 

write! 11,*) 'Electric  Field’ 

write(ll,200)(el(i),el(i+l),el(i+2),el(i+3),el(i+4), 
Xi =0 , 5 , 5  ) 

write!  11, 3  0  0 ) e 1 ( m ) 
write!  11,*)’  ’ 

wr i t e (  1 1 , * )  ’  ’ 

00100  f ormat (  ’ 

’  ,  ’  h  =  ’ , x , e 4 . 2 , x , '  k  =  ' , x , e 6 . 2 , x , ’times’ , x , e 6 . 2 ) 

00200  format!'  ’  , el  1 . 5 , 4 ( 2x , e 1 1 . 5 )  ) 

00300  format!’  ’,ell.5) 

end 
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